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ABSTRACT 

Context. To gain insight on the mass assembly and place constraints on the star formation history (SFH) of Lyman break galaxies 
(LBGs) it is important to determine properly and accurately their properties. 

Aims. We estimate how nebular emission and different SFHs affect parameter estimation and uncertainties. 

Methods. We present a homogeneous, detailed analysis of the spectral energy distribution (SED) of ~ 1700 LBGs from the GOODS- 
MUSIC catalog with deep multi-wavelength photometry from U band to 8 /um, to determine stellar mass, age, dust attenuation and 
star formation rate. Using our SED fitting tool which takes into account nebular emission, we explore a wide parameter space. We 
also explore a set of different star formation histories. 

Results. Nebular emission is found to significantly affect the determination of the physical parameters for the majority of LBGs at 
redshift z ~ 3-6. We identify two populations of galaxies by determining the importance of the contribution of emission lines. We find 
that ~ 65% of LBGs show detectable signs of emission lines, whereas ~ 35 % show weak or no emission lines. This distribution is 
found over the entire redshift range. We interpret these groups as actively star forming and more quiescent LBGs respectively. We find 
that models with constant star formation cannot reproduce the entire range of observed colors, whereas models with nebular emission 
and variable (declining or rising) star formation histories succeed. Other arguments favoring episodic star formation and relatively 
short star formation timescales are also discussed. Taking into account nebular emission generally leads, for a given SFH, to a younger 
age, lower stellar mass, higher dust attenuation, and higher star formation rate, although with increased uncertainties. We find a trend 
of increasing dust attenuation with galaxy mass, and a large scatter in the SFR-M* relation. Our analysis yields a trend of increasing 
specific star formation rate with redshift, as predicted by recent galaxy evolution models. 

Conclusions. The physical parameters of approximately two thirds of high redshift galaxies are significantly modified when nebular 
emission is taken into account. SED models including nebular emission shed new light on the properties of LBGs with numerous 
important implications. 

Key words, galaxies: starburst - galaxies: high redshift - galaxies: evolution - galaxies: star formation 



1. Introduction 

The study of star formation history at high redshift is an im- 



portant issue, with both intensive observationnal (e .g. |Bouwens 
|et al.|2009l [Stark et al.|2009l|Bouwens et al.|2011) and theoret- 
ical work (e.g.|Finlator et al.|2007[|Bouche et al.|2010[|Finlator 



|et al.|2011[[Weinmann et al.|2011| ). While theoretical studies at- 
tempt to reproduce observed trends, as UV luminosity function, 
the lack of available spectroscopic data at z > 3 leads to rely on 
parameter estimation based on indirect indicators as UV J3 slope 
for dust attenuation or UV luminosity for star formation rate (e.g. 
B ouwen s~et al.||2011b or using SED fitting (e.g. |Shapley et al.| 



200 l ||Papovichetal.|2001[ |Verma et al.|2007j |Stark et al.|2 009 ; 



Lee et al.||2011| ). The ability of SED fitting to estimate physi- 
cal parameters has been considerably improved with large multi- 
wavelength surveys, thanks to the commissioning of new space 
and ground-based facilities in the last decades as the Hubble 
Space Telescope, Spitzer or the Very Large Telescope. Indeed, 
an accurate determination of physical parameters is necessary to 
put strong constraints on stellar mass assembly and evolution of 
physical parameters. 

Although nebular emission (i.e. emission lines and nebular 
continuous emission from HII regions) is ubiquitous in regions 



of massive star formation, strong or dominant in optical spectra 
of nearby star forming galaxies, and present in numerous types 
of galaxies, its impact on the determination of physical parame- 
ters of galaxies, in particular at hig h redshift, has been neglected 
until recently (cf. overview in Schaerer & de Barros 2011). 
Several spectral models of g alaxies have indeed included neb- 
ular emission in the past ( e.g. Chariot & Longhetti 200 l\ Fioc & 



Rocca- Volmerang ej 1 997 Anders & Fritze-v. Alvensleben 2003 ; 
Zackriss on et al.||2008| ); however they have not been applied to 
the analysis of distant galaxies. 



Zackriss on et al.| ( |2008| ) showed that nebular emission 
can significantly affect broad-band photometry. Schaerer & de 
Barros (2009 ) included for the first time nebular emission to fit 
SED of a sample of Lyman break galaxies at z ~ 6, and showed 
that nebular lines strongly affect age estimation, since some lines 
can mimick a B aimer break. Ages are strongly reduced, which 
can lead to reconsider star formation rate estimations from UV 
luminosity, since the standard relation used to convert UV lu- 
minosity into SFR is based on the assumption of a constant star 
formation activity during 100 Myr ( |Kennicutt| 1998 ; Madau et al. 
1998 ). The analysis of a sample of z ~ 6-8 LBGs observed with 
HST and Spitzer further demonstrated the potential impact of 



1 



S. de Barros et al.: Impact of nebular emission at high redshift 



nebular emission on t he physical parameters derive d from SED 
fits of high-z galaxies (Schaerer & de Barros 



It has now become clear ([Schaerer & de Barros|20 09, 2010; 



2010). 



|Ono et al.|20T0l|Lidman et al.|2012| ) that nebular emission (both 
lines and continuum emission) must be taken into account for 
the interpretation of photometric measurements of the SEDs of 
star-forming galaxies such as Lyman-alpha emitters and Lyman 
break galaxies - the dominant galaxy populations at high-z. 
Furthermore, as testified by the presence of Lya emission, a 
large and growing fraction of the currently know population 
of star-forming galaxies at high redshift shows emission lines 
([Ouchi etal. |2 008; Star k et al.|2010[|Schaerer et al.|201 lUCurtis^ 
|Lake et al.|2012t|Schenker et al.|2Q12] T 

In parallel a lot of diverse evidence for galaxies with strong 
emission lines and/or strong contributions of nebular emission to 
broad-band filters has been found at different redshifts, e.g. by 
Shi m et al.|(|2011|);IMcLinden et al.] (|201 |Atek et al.| ( [201 1[ ); 



Trump et al.| ( [2011 1 ); | van der Wei et al.| ( |2011| ). It is therefore 



obvious that a systematic study of the impact of nebular emis- 
sion on the physical parameters of LBGs at different redshift is 
urgently needed. This is the main objective of the present work. 

Several other studies included the effect of nebular emission 
(e.g. IPno et al.||2010[ |Schaerer & de Barros||2010[ |Acquaviva| 



etal. 



2011 McLure et al.|2011| ), providing further evidences of 



the impact of nebular emission on physical parameter estima 
tion, leading to lower stellar masses and higher SFRs, possibly 
reconciling expected evolution of increasing specific star forma- 
tion rate (sSFR) with redshift from hydrodynamical simulations, 
with sSFR determined from observations (Bou che et al.|2010| ). 
Indeed, studies which do not consider nebular emission found 
no evolution of sSFR with redshift above z ~ 2 (e.g. |Stark et al. 
|2009l [Gonzalez et al.|2010| ). 

Another issue is that UV p slope seems to indicate that galax- 



ies at high redshift are essentially dust free ( Bouwens et al. 
201 1), while other studies provide evidences for some high red- 



shift obscured galaxies (Verma et al.||2007[ |Yabe et al.||2009[ 
Schaerer & de Barros 2010| ). Furthermore, there is a possi- 



ble trend of increasing dust attenuation with the stellar mass 
( [Schaerer & de Barros|2010|), tre nd already established at lower 
red shift ([Br inchman net al.|2004| ) . 

|Stark et al.| ( |2009| ) was the first attempt of constraining star 
formation history by studying evolution of LBGs sample uni- 
formly selected among different bins of redshift. In this work, 
we use a similar approach, with a large sample of LBGs cov- 
ering four bins of redshift between z ~ 3 to z ~ 6, using an 
up-to-date photometric redshift and SED-fitting tool that treats 
the effects of nebular emission on the SEDs of galaxies. This ho- 
mogeneous analysis provides the main physical parameters, as 
star formation rate, stellar mass, age and reddening. 

Our paper is structured as follows. The observational data 
are described in Sect.[2] and the method used for SED modelling 
is described in Sect. [3[ The results are presented in Sect. [4] and 
discussed in Sect. [5] Section0summarises our main conclusions. 
We adopt a A-CDM cosmological model with Ho=10 km s" 1 
Mpc -1 , £l m =03 and Qa=0.7. All magnitudes are expressed in 
the AB system ( |Oke & Gunn|1983| ) 



2. Data 

2.1. The GOODS Fields 

We focus our analysis on the data from the Great Observatories 
Origins Deep Survey (GOODS). Deta iled descriptions of the 
datasets are available in the literature (Giavali sco et al.|[2004[ 



Santini et a l. 2009 ), so we only provide a brief summary here. 
The GOODS-S and GOODS-N survey areas each cover roughly 
160 arcmin 2 and are centered on the Chandra Deep Field South 
(CDF-S; [Giacconi et al.|2 002) and the Hubble Deep Field North 
(HDF-N; Willi ams et al.|1996 ). Extensive multiwavelength ob- 
servations have been conducted in each of these fields. In this 
paper, we utilize optical imaging from the Advanced Camera 
for Surveys (ACS) onboard the Hubble Space Telescope (HST). 
Observations with ACS were conducted in F435W, F606W, 
F775W, and F850LP (herea fter B435, V606, _ft 7S , 4 50 ) toward 
GOODS-S and GOODS-N ( [Giavalisco et"aLl |2004). The aver- 
age 5cr limiting magnitudes in the v2 GOODS ACS data (07 35 
diameter photometric aperture) are B435 =29.04, V606 =29.52, 
/^ 75 =29.19, and z' 850 =28.54. We also make use of U and R- 
band observations of GOODS-S taken with the ESO V ery Large 
Telescope (VLT) with the VIMOS wide field imager ( |Le Fevre] 
|et al |2003| ), provided by |Nonino et al.| ( |2009| ), with a l<x limiting 
magnitude (17 radius aperture) reaching U ^29.8. 

In the near-infrared, we utilize publicly available deep J, H 
and K-band observations of GOODS-S (PI: C. Cesarsky) using 
the ISAAC camera on the VLT. The sensitivities vary across 
the field depending on the effective integration time and see- 
ing FWHM. Average 5<x magnitude limits (corrected for the 
amount of flux that falls outside of the 17 diameter aperture) 
are / * 25.2, H * 24.7 and K, * 24.7. 

Deep Spitzer imaging is available toward both GOODS fields 
with the Infrared Array Camera (IRAC) as part of the "Super 
Deep" Legacy program (Dickinson et al. in prep). Detail s of 
the observations have been described in detai l elsewhere (Eyles 
|et al.|2005[|Yan et al.|2005||Stark et al.|2007| ) so we do not dis- 
cuss them further here. The 5cr limiting magnitudes of the IRAC 
imaging are ^ 26.3 at 3.6/mi and ^ 25.9 at 4.5/mi using 274 
diameter apertures and applying an aperture correction. 

Optical, near an d mid-infrared photometry are described in 
Sant ini et"aL] ( |2009| ). Non-detections are included in the SED fit 
with Hyperz by setting the flux in the corresponding filter to 
zero, and the error to the l<x upper limit. 



2.2. Dropout selection 

Galaxies at z - 3, 4, 5, and 6 are selected via the presence of the 
Lyman-break as it is redshifted through the U, B435, V606, and 
i' 115 bandpasses, respectively. Selection of Lyman break galaxies 
at these redshifts has now become routine ([Stan way et al. 12003 1 
Gia valisco et al.|2004t |Bunker et al.|2004t |Beckwith et al.|2006[ 
Bouwens et al.|2007] Tln order to ensure a consistent comparison 
of our samples to these previous samples, we adopt color crite- 
ria which are similar to those us ed in |Beckwith et al. ( |2006| ), 
and very similar to those used by Bouwens et al. ( 2007 ). These 
criteria have been developed to select galaxies in the chosen red- 
shift interval while minimizing contamination from red galaxies 
likely to be at low redshift. 

As mentioned in the previous section, galaxies are selected 
from the GOODS-MUSIC catalog v2 ( [Santini et aLl|2009] ). It 
contains 14 999 objects selected in either the zsso band or the K s 
band or at 4.5/mi. 

The selection dropout criteria used in this work are identical 
to those from |Nonino et al.| (|2009) and |Stark etaLl ( |2009| ), with 
an additionnal constraint for U dropout (S/N(U) < 2). This lat- 
ter criteria helps greatly in removing contaminants, but it also 
removes galaxies at the low redshift tail of the U -drop selection. 
This selection leave us 440 U drops, 859 B drops, 277 V drops 
and 66 i drops. 



2 



S. de Barros et al.: Impact of nebular emission at high redshift 



3. Method 

3.1. SED fitting tool 

We use a recent , modified version of th e Hyperz photometric 

( 2000| ), taking into account 



redshift code of Bolzonella et al 



nebular emission (lin es and continua). We co nsider a large set 
of spectral templates (Bruzual & Chariot 2003 ), covering differ- 
ent metallicities and a wide range of star formation (SF) histories 
(exponentially decreasing, constant and rising SF), and we add 
the effects of nebular emission following our method presented 
in |Schaerer & de Barros| ( |2009||2010| ). We account for attenua- 
tion from the intergalactic and the interstellar medium and vary- 
ing redshift. With these assumptions we fit the observed SEDs 
by straightforward least-square minimization. 

In practice we adopt a spectral templates computed for a 
Salpeter IMF ( jSalpeter|1955] ) from 0.1 to 100 M Q , and we prop- 
erly treat the returned ISM mass from stars. Nebular emission 
from continuum processes and lines is added to the spectra pre- 
dicted fro m the GALAXEV models as described in Schaerer & 
de Barros (2009), proportionally to the Lyman continuum pho- 
ton production. The relative line intensities of He and metals are 
taken from |Anders & Fritze-v. Alvensleben (2003), including 
galaxies grouped in three metallicity intervals covering ~ 1/50- 
1 Z Q . Hydrogen lines from the Lyman to the Brackett series are 
included with relative intensities giv en by case B. For ga lactic 
attenuation we use the Calzetti law ([Calzetti et al.||2000j ). The 
IGM is treated following |Madau| ( | 1 995 1 ) . 

To examine the effects of different star formation histories 
(SFH) and for comparison with other studies we define three 
sets of models: 

- Reference model (REF): constant star formation rate, mini- 
mum age t > 50 Myr, and solar metallicity. 

- Decreasing model (DEC): exponentially declining star for- 
mation histories (SFR oc exp(-f/r)) with variable timescales 
r. Metallicity and r are free parameters. 

- Rising model (RIS): rising star formation rate. We use the 
mean rising star-formation history from the simulations of 
Finl ator et al.| ( [201 1] their Fig. 1). To describe this case, we 



assume that SF starts at age=0 and grows by 2.5 dex during 
0.8 Gyr following their functional dependence. After this pe- 
riod we set SFR= 0. Metallicity is a free parameter. 

Furthermore, we define three options concerning the treat- 
ment of nebular emission: 

- No nebular emission. 

- +NEB: including nebular continuum emission and lines ex- 
cept Lya, since this line may be attenuated by radiation trans- 
fer processes inside the galaxy or by the intervening inter- 
galactic medium. 

- +NEB+Lya: including nebular emission (all lines and con- 
tinuum processes). 

SEDs fits to all galaxies have been computed for each of the 
above model sets and nebular options, i.e. for nine different com- 
binations. This allows us to examine in detail the impact these as- 
sumptions/options have on the derived physical parameters, and 
to compare also their fit quality. For the B-drop sample we have 
also tested the effect of the SMC extinction law of Pre vot et aLl 
([1984]). 

In general, the free parameters of the SED fits are: the red- 
shift z, metallicity Z (of stars and gas), the star formation his- 
tory parametrised by r, the age t defined since the onset of star- 
formation, and the attenuation Ay. For the reference model set 



the SFH and metallicity are fixed and the age limited to a mini- 
mum. For the RIS model set, the SFH is also fixed. In all cases, 
we consider z e [0, 10] in steps of 0.1, Ay = 0-4 mag in steps 
of 0.1, and 51 age steps from to the age of the Universe (see 
Bolzonella et al. 2000| ), the SFH of the decreasing models is 
sampled with r = (10, 30, 50, 70, 100, 300, 500, 700, 1000, 
3000, oo) Myr. 

For all the above combinations we compute the x 2 and the 
scaling factor of the template, which provides information about 
the SFR and stellar mass (M*), from the fit to the observed SED. 
Minimisation over the entire parameter space yields the best-fit 
parameters. To determine both confidence intervals (68%) and 
medians for all the parameters, we ran 1000 Monte Carlo simula- 
tions for each object by perturbing the input broadband photom- 
etry assuming the photometric uncertainties are Gaussian. This 
procedure yields the one or two dimensions probability distribu- 
tion function of the physical parameters of interest both for each 
source and for the ensemble of sources. 




Fig. 1. Comparison between photometric redshift and spectro- 
scopic redshifts with in green, black, red and blue the U, B, V, 
/-dropout. Desagreements within 68% confidence limit for U, 
B, V and /-dropout affect 9 objects with very good/good spectro- 
scopic redshift, 6 uncertain and 5 unreliable. Large error bars are 
due to sources with maxima probability at low and high redshift. 



3.2. Redshift selection 

We have compared our photometric redshifts against ob- 
jects with known spectroscopic redshift taken from litterature 



([Vanzella et al. 12005^2 006, 2008; Mign oli et al.|2005|[Szo koly 



et al. 



2004 



2004 



Le Fevre et al.|2004[|Dohertyetal. 2005 ; Wolf etaL 



Daddi et al. 2005 



jCristiani et al.||2 000; Stro lger et al 
20041 |van der Wei et al.]|2005t JRoche et al.||2006l"|Ra vikumar 
et al. 2007; Teplit z et al.|2007| |Xu et al.||2007| |Gronwall et al 



2007; Nilsson et al. 2007; Hathi et al. 2008; Finkelstein et al 



2008) | Yang etaL[2 008 ; Po pesso et al.|2009| ). The result is shown 
in Figure [1J respectively for the U, B, V and i -dropout samples, 
for which 42, 72, 50 and 14 spectroscopic redshifts are available. 



3 



S. de Barros et al.: Impact of nebular emission at high redshift 



To estimate our photometric redshift performance, we compute 
the median Az/(1 + z spe c), with Az, the difference between the 
median z p hot and z spe c- For each sample and model, we obtain 
values from -0.02 to 0.09, with no significant differences among 
models. This result show that we recover redshift with a good 
accuracy, consistent with typical values fo und in other studies 
HWuyts et al.|2009||Hildebrandt et al.|20T0) . 

Combining results of DEC, DEC+NEB +Ly a and 
DEC+NEB models, we find 20 objects among all sam- 
ples, whose median photometric redshift is inconsistent with 
the spectroscopic redshift within the 68% confidence limit. 
The GOODS MUSIC catalog provides quality flags for the 
spectroscopic redshift; among the objects with inconsistent 
redshifts 6 are very good, 3 good, 6 uncertain and 5 unreliable 
spectroscopic redshifts, leading to an estimated 5-11% of 
outliers in our samples. We also obtain objects with large error 
bars, which is due to a double peaked redshift probability 
distribution function with maxima at low and high redshift. 

To eliminate low redshift contaminants and to have the most 
reliable sample at each redshift we proceed with a conservative 
cut: for U, B, V and /-dropouts, we take a lower limit for the 
median photometric redshift of z > 2, z > 3, z > 4 and z > 
5 respectively, as derived from for the DEC, DEC+NEB +Lya 
and DEC+NEB models. We obtain respectively 389 (~ 88%), 
705 (~ 82%), 199 (~ 72%) and 60 (~ 91%) objects. Similar 
criteria applied with REF models leads to larger sample (5 to 
13%) and with RIS models to similar samples with a maximum 
variation of the size of the sample of 4%. We notice that with 
this final selection, there is a significant overlap between the U 
and B dropout with 96 objects in both samples. 

Table Q] shows median redshifts and 68% confidence lim- 
its for REF, DEC and RIS model. Considering +NEB+Lya and 
+NEB models induces median redshift variation not higher than 
0.1. 



Table 1. Median redshift values and 68% confidence limits of 
the final samples. 



REF 



DEC 



RIS 



t/-dropout 3.32^ 
5-dropout 3.88^J 



V-dropout 
/-dropout 



4.94" 
6.00" 



M 

-0.29 



o oa+0.26 

-Q-14 
o 7Q+0.44 

J - ' -0.39 

4.81+ 059 

6.O0M 



3.30 + ^5 



3.78 
4.81! 
6.00 + 



-All 



4. Results 

4.1. Two LBG categories revealed 

As this is the first time large samples of LBGs are analysed with 
SED fits including the effects of nebular emission, we have ex- 
amined if this leads to better fits and by how much. Our main 
result from this comparison is that typically 60-70% of the 
galaxies are better fit with nebular emission (option +NEB or 
+NEB+Lya) than without. This is found quite independently 
of the adopted SFH and for all samples, i.e. from z ~ 3 to 6. 
In other words -35% of the objects are best fitted without tak- 
ing account of nebular emission. This fraction is independent 
of properties such as the absolute UV magnitude M1500 or the 
number of filters available. Furthermore, all SF histories (REF, 
DEC, RIS model sets) yield approximately the same percent- 
ages (30%-39%). Finally, this is not only a statistical property, 




-L5 -L0 



3,6/jm-4.5/jm [mag] 



Fig. 2. 3.6yum-4.5yum color histogram for a sub-sample of z 
e [3.8, 5] objects. In blue, best fitted objects with nebular emis- 
sion and red, best fitted objects without nebular emission, both 
for decreasing SF. 



but the vast majority of objects can be assigned to such a cate- 
gory. Indeed, among the U, B, V and /-dropouts, we find 68%, 
71%, 71% and 77% common objects best fitted without nebu- 
lar emission and 85%, 80%, 94% and 88% common objects best 
fitted taking account nebular emission. 

Since this distinction in two groups is fairly model- and 
redshift- independent there must be a physical explanation for 
it. The easiest and most natural explanation is found when con- 
sidering a sub-sample of objects over a restricted redshift inter- 
val. Indeed, since Ha is a strong line at 656.4 nm (rest-frame) 
and few strong lines are found longward of it, this line must af- 
fect the 3.6-4.5 jim color for objects between z=3.8 and z=5 (cf. 
Shim et al.|2011| ). We therefore selected ^-dropout objects with 
available 3.6/mi and 4.5/mi data (excluding non-detections) and 
with median redshift between 3.8 and 5. We obtain a subsample 
of 303 objects, for which again -35% are better fit without neb- 
ular emission. This should thus be a representative subsample 
of all galaxies studied here. Figure [2] shows that the objects best 
fitted with nebular emission have a systematically bluer 3.6/mi- 
4.5yum color than those better fit excluding nebular effects. This 
shows that objects better fit with models accounting for nebular 
emission do indeed show strong Ha emission lines. This is not 
a trivial finding, since these models also allow for ages/SF his- 
tories where nebular emission is absent/insignificant. We there- 
fore conclude that the objects best fit with models accounting 
for nebular emission (~ 60-70%) correspond to galaxies with 
"strong" emission lines, whereas the rest shows few or no dis- 
cernible signs of emission lines. 

It is inter esting to note that both Yabe et al. (2009) and Shim 

eTaL] pOTT]) also found ~ 70% of a sample of LBGs at z ~ 5 



and 4 respectively with a 3.6 fim excess, which they explained 
by the presence of Ha emission. While their samples include ~ 
100 (70) such galaxies, we have - 300 galaxies at z ~ 3.8-5 
with direct empirical evidence for strong Ha emission. In ad- 
dition our SED fits suggest that a similar percentage of objects 
with "strong" emission lines exists over the entire examined red- 
shift range (z ~ 3-6). Below we will show that this interpretation 
is perfectly consistent with the differences found for the physi- 
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cal parameters of these galaxies, and we will propose a physical 
explanation for the existence of these two groups of LBGs. 
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Fig. 3. Comparison between DEC model x 2 anc * 
DEC+NEB +Lya/+NEB model xl- R ed triangle: DEC model 
best fit and blue square: DEC+NEB +Lya/+NEB best fit. Blue 
squares and red triangles underly the two LBGs categories, 
respectively "strong" nebular emitters and "weak" nebular 
emitters. 



4.2. Fit quality and constraints on the star formation histories 

To compare the fit quality of the different models, we compare 
values of x 2 r (X 2 value divided by the number of filters minus 
1). At each redshift SEDs are systematically better fitted with 
RIS and DEC model sets (considering or not nebular emission) 
in comparison with REF model sets. x 2 r values are on average 
20-40% lower for RIS model sets and 25-40% lower for DEC 
model sets, which also shows that DEC model sets fit slightly 
better than RIS model sets. As shown in Figure [3] for declining 
SF and all redshifts, "strong" nebular emitters show a large im- 
provement in x 2 when they are fitted with models considering 
nebular emission (x 2 ~30 to 55% lower on average). At the op- 
posite, "weak" nebular emitters show a slight improvement of 
their x 2 when they are fitted without nebular emission (10 to 
30% on average). Models with nebular emission are able to pro- 
vide significantly better fits for "strong" nebular emitters and a 
more or less similar fit quality for "weak" nebular emitters. 

For the subsample of objects with redshifts z ~ 3.8-5 dis- 
cussed above, the observed 3.6yum-4.5yum provides also an inter- 
esting constraint on star formation history and nebular emission. 
Indeed, as shown in Fig. |4| models without nebular emission 
are unable to reproduce the range of observed 3.6yum-4.5yum col- 
ors at z ~ 4. The REF+NEB model (constant star formation) is 
also unable to reproduce observations while the DEC+NEB and 
RIS+NEB models provide fits fully consistent with this color 
(+NEB+Lya option leads to similar results), with DEC+NEB 
providing even better results than RIS+NEB. Abandoning the 
age limitation (age > 50 Myr) for the REF+NEB model 




-0.5 0.5 1 
3.6^m-4.5/xm obs 



0.5 0.5 1 

3.6^m-4.5/im obs 



Fig. 4. Comparison between observed 3.6yum-4.5yum color and 
best fit color for a sub-sample of z 6 [3.8, 5] objects for SFHs 
with and without nebular emission. 



would allow this model to provide fits consistent with observa- 
tions. We need to have very large Ha equivalent widths (EWs) 
to reproduce the color of "strong" nebular emitters. This can 
be obtained for any SFH with young ages (median of ~ 20 
Myr for DEC+NEB and RIS+NEB). On the other hand for 
"weak" nebular emitters small EW(Ha) are necessary to repro- 
duce their red color. Since both for rising and constant SFs no 
physical process is able to suppress effectively nebular emis- 
sion, DEC+NEB (+Ly a) is the model which best fits both LBG 
categories. However, considering errors in color estimation, 
RIS+NEB(Lya) and REF+NEB(+Lya) provide also acceptable 
fits for "weak" nebular emitters. Further constraints and results 
on the timescales of the exponentially declining star formation 
histories are discussed in Sect. I4.3.5I 

Finally, among objects better fitted with nebular emission we 
have found a shift between those best fit with or without Lya 
(i.e. between +NEB and +NEB+Lya option) with redshift, in 
the sense that higher-z galaxies favor a larger fraction of objects 
with Lya emission. This shows that SED fitting is also sensitive 
to Lya emi ssion, a finding we ha ve demonstrated and discussed 
in detail in Schaer er" "t al.| ( |20lT] ). 

4.3. Physical properties 

We now turn to discuss the main physical properties (stellar 
mass, SFR, age, attenuation, plus the star formation timescale 
where appropriate) of the LBGs and their dependence on model 
assumptions. The median values and uncertainties of the phys- 
ical parameters derived in several bins of UV magnitude Mi 500 
for all our samples, and using all nine combinations of model 

and IA.3I For each 



A.2 



assumptions, are listed in Tables |A.l 
physical parameter we will now describe the median properties 
and their model dependence, explain their origin, and compare 
the behaviour of the individual values. Furthermore we will ex- 
amine possible correlations between derived and observed pa- 
rameters. To do this we here chose the largest subsample, con- 
sisting of 705 B-drop (z ~ 4) galaxies, since overall the same 
trends/differences are found at all redshifts (except stated other- 
wise). In Sect. |4.4| we then discuss the redshift evolution of the 
physical parameters and their model dependence. 
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4.3.1. Absolute UV magnitude 

In what follows the absolute UV magnitude Mi 500 refers to the 
absolute magnitude at 1500 A. To determine it for each object, 
we use the integrated SED flux in an artificial filter of 200 A 
width centered on 1500 A. Using the V-band magnitude for U- 
drop, /-band for #-drop, and z-band for V and /-drop samples and 
spectro scopic redshift wh en availaible, to estimate the UV mag- 
nitude ( [Stark et aL]|2009| ) leads to no significant difference on 
M1500, except for one Z?-drop, two V-drop and one /-drop galaxy. 
These are in fact objects with spectroscopic redshift identifica- 
tion at low redshift which pass our selection, and as they repre- 
sent less than 1 % (slightly more for /-drop) of each sample, we 
consider that they can not alter significantly our conclusions. In 
passing we note that the number of objects in each UV magni- 

[|A.3|can somewhat change 



tude bin listed in Tables A. 1 A.2|and 



from one model to another, mostly due to small differences in 
photometric redshifts. 

4.3.2. Age 

Overall, the ages of individual z ~ 4 LBGs derived from the 
different models span a wide range, typically from ~ 4 Myr (if 
no lower limit is specified) to ~ 1 Gyr, the maximum age at this 
redshift (see Fig. [5]). A wide age range is found for all 9 nine 
model sets. 

The individual ages and the resulting median age of the sam- 
ple depend qui te strongl y on the model assumptions. As can be 
seen in Tables A.l - A.3 the assumption of declining (DEC) or 
rising (RIS) star formation histories lead - for models without 
nebular emission - to median ages younger by a factor 5-10 
compared to constant star formation (REF). These differences 
are much larger than the typical age uncertainty (~ 0.15 dex) 
found for the REF model. The reason is that with a constant SF, 
galaxies keep a high UV rest- frame flux and if the observed SED 
shows the presence of a B aimer break, older ages are necessary 
to have a sufficient population of evolved stars to reproduce the 
observed break. For a declining or rising SF, much younger ages 
are sufficient to obtain a suitable evolved population. Indeed the 
observed median ratio of the optical/UV flux i s high and con- 
sistent with a Balmer break at each redshift (see jGonzalez et al.] 
|20T0l|Lee et al.|2011| ). 

As already shown in Schaere r & de Barros ( [2009 2010 ), 
emission lines (mostly [O III] /U4959, 5007 plus HJ3) can mimic 
a Balmer break for a constant or decreasing SF, and our results at 
Z ~ 3-5 demonstrate this effect for a large sample: the inclusion 
of nebular emission (with or without hya) leads to systemati- 
cally younger median ages for the REF and DEC model sets (e.g. 
by a factor 2-5 for the REF model; see Fig. [5]). This is not the 
case for the rising star formation history: RIS+NEB+Lya and 
RIS+NEB models lead to older or similar ages, for reasons ex- 
plained below. Furthermore, as shown in Schaerer & de Barros 
( 2010| ), taking into account nebular emission increases the un- 
certainty on the age, especially for REF and RIS model sets, 
from 0.15 to 0.27 dex and from 0.39 to 0.55 dex respectively 
(~ 0.4 dex for DEC model sets). However, these uncertainties 
are sufficiently small to distinguish differences in age among the 
three SFHs (constant, decreasing and rising) and between mod- 
els with/without nebular emission. 

To discuss the effect of nebular emission on age it is useful 
to consider the two LBG categories, "weak" and "strong" line 
emitters, separately. Median ages of the "weak" nebular emit- 
ters increase systematically when we take into account nebular 
emission, for all considered star formation histories. This is nat- 




7 8 9 7 8 9 

Log Age REF [yr] Log Age REF [yr] 



Fig. 5. Composite probability distribution of age for REF model 
and age for all other models at z ~ 4. The points overlaid show 
the median value properties for each galaxy in the sample, black 
dots for "weak" nebular emitters and white dots for "strong" 
nebular emitters. The overlaid contour indicates the 68% inte- 
grated probabilities on the ensemble properties measured from 
the centroid of the distribution. 



ural, since in all cases the equivalent widths of the optical emis- 
sion lines decrease with time. Hence increasing ages minimizes 
their contribution. For the RIS model, the ages of "weak" neb- 
ular emitters are much more increased than for REF and DEC 
models, by a factor 10-20 at z ~ 4 (compared to 2-5 times for 
REF and DEC models). This strong increase outweights the age 
decrease of the "strong" emitters, which explains why the me- 
dian age of the whole population increases or does not change, 
when we take into account nebular emission for the rising star 
formation history. 

Quite generally, and for all SFHs, we find for "strong" neb- 
ular emitters median ages, which are systematically younger 
when nebular emiss ion is included, because of a stron ger effect 
Schaerer & de Barros| ( [2009l[20T0l ) and Con- 
ines mimick a Balmer break. Only at z ~ 6 



already discussed in 
firmed here: strong 



and for DEC+NEB and RIS+NEB models do we find older ages. 
Balmer break seems to be absent in these galaxies, which are fit 
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with a strong UV flux (models with no nebular emission) and/or 
strong Lya emission (+NEB+Lyor). For DEC or RIS models, 
this implies a young population, and for DEC/RIS+NEB+Lya 
too. For DEC/RIS+NEB models, the lack of Lya emission must 
be balanced by an increase of UV flux, which implies a younger 
stellar population, but a younger population leads to strong lines 
which mimick a strong B aimer break, this is why age cannot be 
as young than for +NEB+Lya models. 

In absolute terms, the derived ages of "weak" and "strong" 
nebular emitters are fairly similar when constant SF is assumed. 
For decreasing and rising SF, the median ages of "strong" neb- 
ular emitters are systematically older than those of the "weak" 
emitters for models without nebular emission, while the oppo- 
site is found for models with nebular emission. In the former 
case (no nebular emission) the "strong" emitters are fitted with a 
strong B aimer break, in the latter case with strong lines. In any 
case, the age differences seem to confirm an intrinsic difference 
between these two LBG categories. 




4.3.3. Reddening 

The attenuation of individual z ~ 4 LBGs derived from the SED 
fits ranges from Ay = to a maximum of ~ 1.5-2 mag for few 
objects, as shown in Fig. [6] Although all model sets yield a sim- 
ilar range of attenuation, relatively large systematic differences, 
which we will just discuss, are found between them. 

Considering or not nebular emission for a given SFH lead 
to variation no larger than 0.2 mag. Comparing constant SF with 
other star formation histories, the median Ay is higher for declin- 
ing SF (+0.2 mag) and for rising SF (+0.5 mag), which can be 
partially explained by the well-known degeneracy between red- 
dening and age, since declining and rising SF lead to younger 
median ages than constant SF (cf. above). Furthermore, since 
young stars always dominate the UV flux for rising star forma- 
tion histories, a higher attenuation is required to fit the observa- 
tions, c ompared to models with c onstant SF, as already pointed 
out by Scha erer & Pello| ( |2005| ). Taking into account nebular 
emission leads on average to variations no larger than 0. 1 mag 
in the reddening for constant SF; a slight increase is found when 
the Lya line is included (DEC+NEB+Lya model, assuming the 
maximum case B intensity for Lya), which is explained by the 
additional flux from the Lya line. 

On the other hand, for rising star formation rate, taking into 
account nebular emission leads to a systematically lower median 
Ay, since the contribution of strong nebular lines in the opti- 
cal and infrared (rest- frame) leads to redder SEDs. These results 
have to be taken with some caution since typical errors on Ay 
for individual objects are 0.1 mag for constant SF and ~ 0.2 
mag for decreasing and rising SF, which does not allow a clear 
distinction for example between REF and DEC models. 

The inclusion of nebular emission with constant star forma- 
tion history does not lead to any significant change in dust red- 
dening for "weak" nebular emitters. On the other hand, for de- 
clining and rising SF, the extinction decreases strongly (by ~ 
-0.2 to -0.5 mag in Ay) when we consider nebular emission. 
Since these objects seems to have intrinsically no discernible 
signs of emission lines, models including nebular emission fit 
by minimizing equivalent widths, and so, by minimizing SFR 
and UV flux. Model set based on constant star formation history 
(REF/REF+NEB/REF+NEB -\-Lya models) do not allow suffi- 
cient variations to modify the reddening estimation. Median Ay 
between z ~ 3 and ~ 5 for "strong" nebular emitters increases 
when we consider decreasing star formation (+0.2 to +0.5 mag), 
while it decreases for rising SF (-0.1 to -0.2 mag). Indeed, fitting 
"strong" nebular emitters with DEC model leads to increase the 
purely stellar continuum contribution in the optical/IR and there- 
fore increase the age and decrease SFR and UV flux. Rising SF 
leads to very large equivalent widths (in comparison with other 
SFHs) in optical and near-infrared (rest- frame), which redden 
SEDs and allow lower values of Ay. 

At z ~ 6, effects of modeling with nebular emission are dif- 
ferent: for decreasing and rising model sets, considering nebular 
emission leads to a decrease of median Ay for both "weak" and 

"strong" nebular emitters, respectively with 0.6 and 0.2 

mag. SED fits with nebular emission lead to an important con- 
tribution of nebular lines longward UV, and strong lines being 
associated with strong UV flux, an additionnal amount of dust 
attenuation is required. 



Fig. 6. Same as Figure^ for Ay. 



4.3.4. Reddening and UV-slope 

Since often the observed UV slope J3 is used to measure the atten- 
uation in LBGs, it is interesting to examine how the attenuation 
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Fig. 7. UV-continuum slope ft vs Ay for all models at z~ 4. Each cross shows the median value properties for each galaxy, blue for 
"strong" nebular emitters and red for "weak" nebular emitters. Black squares are median values by bins of 0.1 Ay mag. The red 
dashed line is the relation be tween extinctio n ( Calzetti et al.|2000" ) to a given f3 if the base spectrum is a young star-forming galaxy 
of constant star-formation ( [Bouwens et al.||2009| ). The red line is a linear fit among the whole composite probability distribution 
function. 



derived from the SED fits based on the different model assump- 
tions compares with /3. Such a comparison is shown in Fig. [7] 
for the nine model sets applied to the 5-drop sample. The UV 
slope has been determ ined using the same filters and relations as 
Bouwens et al.| ([2009 ). Figure [7] shows that there is a significant 
trend of increasing J3 with Ay, as expected, albeit with a large 
scatter. We have done linear fits to the 2D composite probability 
distribution function, which yields the mean relations indicated 
in the plot by red lines. For compariso n the "standard " relation 
between J3 and Ay, here taken from B ouwens et al.| ( |2009| ), is 
also shown. As expected, our relations agree well with the "stan- 
dard" one for models assuming constant star formation and ages 
> 50 Myr (REF model sets), since this corresponds to the main 
assumption made to derive the standard /3-reddening relation. 
For a given SF history, the differences between the three op- 
tions with/without nebular emission can be explained by the be- 
haviour of A v discussed above. Since all models with declining 
and rising star formation histories yield somewhat higher red- 
dening on average (cf. above), a relation steeper than the stan- 
dard one is found. Since the obtained relations are fairly similar 



we can combine them to obtain the following mean relations be- 
tween p and Ay for the three cases, ([I]) modeling without nebular 
emission, ^ +NEB+Lya and ([3} +NEB: 



1.54 x OS + 2.54) 
1.47 x OS + 2.49) 
1.09x00 + 2.58) 



(1) 

(2) 
(3) 



The last relation (Eq.|3]) is probably the most appropriate one, 
since it combines the models which best fit the data (i.e. models 
including nebular emission but no strong Lya for the majority of 
the galaxies, cf. |Schaerer et al.|201 1| ). It should be reminded that 
these relations assume a Calzetti attenuation law; to translate this 
into the color excess one has E(B - V) = A v /R v , where R v = 
4.05. In short, from our new relations derived from a subsample 
of 705 B-drop galaxies using various star formation histories we 
find that LBGs with a given UV slope have higher attenuation 
than derived from the commonly used fi-Ay relation. For typical 
average UV slopes of J3 2.2 (- 1 .7) found for faint (bright) z ~ 
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A [A] 

Fig. 8. Theoretical SEDs in the rest- frame, normalized at 2000A, 
for models with EW(Ha) = 500 A (solid lines) and 100 A 
(dashed lines) and different star formation timescales (constant 
SF in red and r = 10 Myr in blue). For the models with 
EW(Ha) = 500 A the ages are 52 Myr and 16 Myr respectively, 
for EW(Ha) = 100A, 2. 1 Gyr and 35 Myr. 



4 galaxies, this translates to an increase of the UV attenuation by 
a factor ~ 3. 



4.3.5. Star formation timescale 

For the models with exponentially declining star formation his- 
tories (DEC models), which are found to provide the best fits for 
the majority of objects - although not always by large margins 
- it is of interest to examine the resulting timescales r and the 
uncertainties on this quantity. As a reminder, the DEC model set 
considers 10 different star formation timescales r € [10,3000] 
Myr plus the limiting case of r = oo corresponding to constant 
star formation. 

For all options with/without nebular emission (DEC, 
DEC+NEB+Lya and DEC+NEB) we find median values of r 
between 10 and 300 Myr in the different UV magnitude bins (cf. 
Table A.2). Models including nebular emission favour on aver- 
age somewhat shorter timescales than those without. 

Although the timescales found are relatively short compared 
to the dynamical timescale - t& yn ~ 40 Myr at z ~ 4 - the uncer- 
tainties on r are quite large. 

What constrains the SF timescales and why are short 
timescales preferred? Since SFRoc exp(-f/r) two (or more) ob- 
servational constraints are need to determine both the age t and 
timescale r. Lets us first consider the case of models including 
nebular emission and examine galaxies with z ~ 3.8-5. In this 
case we find that t and r are mostly constrained by a combination 
of the (3.6-4.5) pm color tracing EW(Ha), and by a UV/optical 
(rest- frame) color. This works as follows. As already discussed 
above, a measurement of the 3.6-4.5 pm color for galaxies at 
Z ~ 3.8-5 reflects the Ha equivalent width. However, it is well 
known that a given equivalent width can be obtained with differ- 
ent values of t/r (cf. Fig. 10). To illustrate this we show in Fig. 
[8] the predicted SEDs for galaxies with the same EW(Ha) = 100 
(500) A but different SF timescales (r = 10 Myr and oo). It is ob- 
vious that the main feature allowing to lift this degeneracy is the 
ratio of the UV/optical flux. At z ~ 4 this ratio is e.g. reflected 




-2 -10 1 2 3 4 
^ssoLF —4.5pm [mag] 



Fig. 9. Z85olp - 4.5/zm observed color histogram at z ~ 4. In red, 
best fitted objects with r > 100 Myr and blue, best fitted objects 
with r = 10 Myr, for the three models with decreasing SF. The 
blue, red and black dashed lines show respectively median color 
for objects best fitted with r = 10 Myr, with r > 100 Myr and for 
the whole sample. For the two subsamples, KS test for the three 
models gives p = 7.8. 10" 8 , p = 2.6. 10" 17 and p = 1.3.10" 16 
respectively for DEC, DEC+NEB +Ly a and DEC+NEB models, 
showing that the two subsamples are not drawn from the same 
population. 



by the zssolp - 4.5/zm color, whose evolution with t and r is also 
shown in Fig. 10 Therefore these two colors will provide good 
constraints on t and r, provided they are not strongly affected 
by reddening (cf. below). Indeed within typical 68% error bars, 
these two extreme SFHs can be discriminated by their color and 
EW(Ha) for t ^ 20-30 Myr. For younger ages, the uncertainties 
in both EW(Ha) and color do not allow a clear separation. A 
posteriori we can verify that the objects best fitted with "long" 
timescales do indeed statistically differ from those with "short" 
timescales. Indeed Fig. [9] shows a statistically significant differ- 
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Fig. 10. Evolution of EW(Ha) (blue) and zssolp - 4.5yum color 
(red) with age for r = 10 Myr (solid lines) and for r = oo (dashed 
lines). Typical error bars, shown on the left for EW(Ha) and on 
the right for the Zssolp ~ 4.5yum color, have been derived from 
error estimation on measured fluxes at z ~ 4. The effect of red- 
ddening (A v = 0.5) on zssolp ~ 4.5/rni is shown with the red 
arrow. 



ence between galaxies best fitted with r = 10 Myr, which are 
bluer in z^solp ~ 4.5//m, and r > 100 Myr showing redder colors. 

We also have to consider dust attenuation, which can in- 



crease Z850LP -4.5/zm as shown in Fig.[10]and introduce a degen- 
eracy between SFHs. However, considering the observed median 
color and dust attenuation for the REF and DEC model sets, we 
would expect to have higher extinction for REF model sets than 
for DEC model sets, but we find the opposite with median Ay 
larger by 0.1-0.2 mag for DEC models. Furthermore, Figure 10 
shows that it is more difficult to reproduce red color with a con- 
stant SF than with a declining SF. This shows that dust attenu- 
ation seems to be only weakly correlated with zzsolp - 4.5yum 
color at z ~ 4, allowing us to conclude that the ratio of UV/near- 
IR flux and equivalent widths of different emission lines (mainly 
Lya, OIII and Ha) drive the choice of r for models with declin- 
ing star formation. Despite the possibility to discriminate two 
extreme SFHs like constant SF and decreasing SF with r = 10 
Myr, large uncertainties on r estimates are found (from 0.7 dex 
for DEC+NEB+Lya model to ~ 2 dex for DEC+NEB model), 
preventing us from deriving strong constraint on r. These uncer- 
tainties come mainly from reddening; indeed fixing reddening to 
an arbitrary value (A v - 0) leads also to a low median r (< 300 
Myr) but with lower typical uncertainties, from ~ 0.4 dex for the 
DEC+NEB +Ly a model to ~ 0.6 dex for the DEC model. 

Interestingly, our models with declining SF histories indicate 
a possible increase of the timescale r with UV luminosity, and 
hence also with stellar mass, but only for "strong" nebular emit- 
ters. It is tempting to suggest that this could be due to a decrease 
of the feedback efficiency with increasing galaxy mass, since 
the star formation timescale is likel y related to the dynam ical 
one and modulated by feedback (e.g. Wyit he & Loeb|20lT] ). We 
do not see any evolution of the timescale for the "weak" nebu- 
lar emitters, which is easily explained by weaker constraints due 
to the absence of strong distinctive features. This explains why 
the trend of r with UV magnitude cannot be seen in Table |A.2| 
where the combined data for the entire sample is listed. 
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Fig. 11. Same as Figure 



4.3.6. Stellar mass 

Stellar mass is generally considered the most reliable parameter 
estimated by SED fitting, since relatively small differences are 
found when varying assumptions like the star formation history 
or dust extinction. Finlator et al.| ( |2007| ) estimates differences 
due to differ ent assumptions on SFHs typically not higher than 
0.3 dex and |Yabe et al. ( 2009| ), adding effects of metallicities 
and extinction laws, estimates differences to be not higher than 
-0.6 dex. Our comparison of stellar masses based on the differ- 



ent model assumption is show in Fig. 1 1 and in Tables A.1-A.3 
Our results confirm earlier results on the dependence of stellar 
mass with the assumed SFH. Indeed, from z ~ 3 to z ~ 6 the me- 
dian stellar masses do not differ by more than -0.3 dex between 
different star formation histories, even with metallicity included 
as a free parameter. With respect to models for constant star 
formation and without nebular emission (REF) all other mod- 
els and options (+NEB or +NEB+Lya) lead systematically to 
lower stellar masses, with differences exceeding the typical un- 
certainties of ~ 0.15 dex found for the REF model. Taking into 
account nebular emission for constant SF leads to differences of 
same order as the uncertainties, not larger than ~ 0.2 dex. For 
models with decreasing or rising SF and nebular emission in- 
cluded, we obtain stellar masses lower by ~ 0.4 dex on average 
compared to the model with constant star formation. 

As could be expected, the masses of galaxies from the cate- 
gory of "strong" nebular emitters vary in the same manner, but 
more strongly than those of the "weak" emitters, when different 
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Fig. 12. Composite p robability distribution of M1500 and M+ for all the models at z ~ 4. The black dashed line represents the M*- 
M1500 trend found by Gonza lez et al. ( 2011 ) and black dotted lines shows a scatter of ±0.5 dex. The solid red line shows a linear fit 
established taking into account the whole composite probability distribution. The points overlaid show the median value properties 
for each object in the sample, black dots for "weak" nebular emitters and white dots for "strong" nebular emitters. The overlaid 
contour indicates the 68% integrated probabilities on the ensemble properties measured from the centroid of the distribution. 



SF histories and options on the treatment of nebular emission are 
taken into account. Typically masses decrease by ~ 0.4-0.9 (0.2) 
dex for "strong" ("weak") nebular emitters between constant SF 
and the other models. With declining and rising SF histories and 
SEDs including nebular emission, the categories of "strong" and 
"weak" nebular emitters show distinct median masses — the for- 
mer being less massive than the latter - with differences up to 
~ 0.6 dex at z ~ 6. In contrast, models with constant SF and 
nebular emission do not yield significantly different masses be- 
tween these LBG categories. This is due to the more restricted 
"dynamical range" allowed by these models (see e.g. Fig.[T0|). 



In Figure 12 we show the stellar mass-Mi5 o relation found 
for all our models at z ~ 4. For constant star formation with 
or without nebular emission we find, as expected, a relation in 



good agreement with that established by Gonza lez et al. ( 2011| ), 
within a scatter of ±0.5 dex. Indee d, our R EF mod el is b ased 
on assumptions similar to those of Gonza lez et al.| ( |2011| ), ex- 
cept for the metallicity (we assume Z = Z Q when they assume 
Z = O.2Z ) and for the minimum age (we assume 50 Myr when 
they assume 10 Myr). Solar metallicity leading to ~ 0.06 dex of 
increase in mass in comparison with 0.2 Z© and a higher mini- 
mal age leading to higher lower bound in M* , the slight offset 
of our stellar mass-M 1500 ^ relation is easily explained. Overall, 
although differences exist in the M* -M\ 500 relation obtained be- 
tween different model sets, our mean relation remains (within 
the scatter of ±0.5 dex) fairly similar to the relation derived by 
Gonzalez et al. ( 2011| ), regardless of our model assumptions. 
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Fig. 13. Composite probability distribution of M+ and age for 
some models at z ~ 4. The points overlaid show the me- 
dian value properties for each object in the sample, black dots 
for "weak" nebular emitters and white dots for "strong" nebu- 
lar emitters. The overlaid contour indicates the 68% integrated 
probabilities on the ensemble properties measured from the cen- 
troid of the distribution. 



A correlation between stellar mass and age is found for all 
the models, as shown in Figure [T3| Since age estimation depends 
of the mass to light ratio, this relation is trivial if we describe star 
formation history with a monotonic function. Indeed, luminos- 
ity is fixed by both measured fluxes and redshift, and NIR data 
putting strong constrains on stellar mass estimation, the stellar 
mass-age relation simply reflects the increase of mass to light 
ratio with age. As previously explained, different assumption on 
the SFH lead to different trends between "weak" and "strong" 
nebular emitters, when analysed with SEDs including nebular 
emission: for variable star formation histories (both rising or de- 
clining) "weak" nebular emitters are found to be older and more 
massive on average than "strong" nebular emitters. In contrast, 
when constant star formation is assumed, the physical properties 
of the two populations do not differ. 
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for Ay and M*. Bottom: rela- 



In Figure 14 (top panel), we show the relation between the 
dust attenuation A v and stellar mass, for a selected model set. |2007[|Maiolino e t al. 2008). 
For all models a similar trend is found with the median Ay in- 
creasing with galaxy mass, and a wide range of attenuations al- 
lowed between zero and -1.5 mag. 



tion between stellar mass and reddening for DEC+NEB model 
at z ~ 4. Blue dots for galaxies with median age < 10 7 years, red 
squares: 10 7 < age < 10 8 , yellow upward triangles: 10 8 < age 
< 10 9 and black downards triangles: age > 10 9 years 



ties in both studies. While this trend could be explained by the 
age-reddening degeneracy, fixing age at a given value lead to no 
change. The most natural explanation of this trend is probably 
that the dust attenuation is related to the stellar mass-metallicity 
relation (cf . |Tremonti et al |2004[|Erb et al.|2006[|Finlator et al] 



4.3.7. Star formation rate 



The bottom panel of Figure [14] allows a better understand- 
ing of the correlation between stellar mass and dust reddening, 
since it shows that the dispersion comes mainly from the age 
scatter. Indeed, for a given range of age, we find a clear trend 
of increasing extinction with increasing stellar mass. This trend 



has a l ready been highlighted at lower redshift (eg. |Buat et al. 
20081 [Burgarella et al,||2007| [Daddi et al.||2007[ |Reddy 



2005 



et al. 



|2006M2008t |Sawicki 



2010| ). Other studies (e.g. Bouwens et a, 



2012]) and it seems a 



so to be ob- 



served at higher redshift ( Yabe et aL[2 009 ; Schaerer & de Barros 



the existence of this trend, while the resu 



2009) also suggest 
from [Dominguez 
et al.| ( [20T2] ), revealing also a similar trend up to redshift 1.5 
is clearly compatible with our result, but with large uncertain- 



The star formation rate (defined here as the instantaneous value 
at the age t) depends quite strongly on the model assumptions, as 
illustrated in Fig. [15] For constant star formation, the inclusion 
of nebular emission leads to somewhat higher SFRs on average, 
due to the younger age, which requires a higher attenuation. The 
largest differences (up to ~ 1 dex) with respect to constant SFR 
models are obtained with exponentially declining star formation 
histories. The reason for such differences is obviously due to the 
variations in the UV output with time and due to the allowance 
of young ages (< 50 Myr), which also imply a higher attenuation 
on average (cf. above). 

An interesting feature of the declining SF histories is that it 
also allows for SFRs lower than those derived using the canoni- 
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Fig. 15. Same as Figure^for SFR. 



cal calibrations assuming constant SFR. Typically, galaxies with 
t/r < 1 , have higher SFR (up to 1 dex), if the timescale is short 
enough to diverge significantly from constant SF. Galaxies with 
t/r > 2 are much more quiescent [] and have lower SFR (up to 
1 dex), and at last, intermediate galaxies have SFRs consistent 
with results from the REF model. This larger "dynamic range" 
may well be physical, as indicated by the existence of "strong" 
and "weak" nebular emitters, as we will discuss below. Models 
with rising star formation histories lead to the highest SFRs, 
since their SED is always dominated by young stars, which im- 
plies a narrower range of UV- to-optical fluxes, hence requiring 
on ave rage a higher atten uation than for other SF histories (cf. 
above, Schae rer & Pel lo 2005 ). For rising SF, with ages a bove 
~ 10 8 yr, all galaxies follow the canonical relation ( [Kennicutt 
|1998| ). Below this age, for the same reason as for decreasing SF, 
the SFR estimated by SED fitting is higher (regardless of dust 
reddening). 

Taking into account nebular emission does not lead to any 
significant change in the median SFR for constant star forma- 
tion, while for decreasing SF the median SFRs are lower (mainly 



1 By quiescent galaxies we here mean objects with lower SFR, not 
galaxies with SFR~ 0. 
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Fig. 16. Composite probability distribution of Mi 50 o and SFR for 
the REF (top), DEC+NEB (centre) and RIS+NEB model (bot- 
tom) for the sample at z ~ 4 as determined for each galaxy from 
our 1000 Monte Carlo simulations. The points overlaid show the 
median value properties for each object in the sample, black dots 
for "weak" nebular emitters and white dots for "strong" nebular 
emitters. The overlaid contour indicate the 68% integrated prob- 
abilities on the ensemble properties measured from the centroid 
of the distribution. The dashed line represents the Kennicutt re- 
lation ( |Kennicutt|1998| ). 



for +NEB+Lya) or equal (mainly for +NEB), and for rising SF 
median SFRs are systematically lower, typically by a factor ~ 2. 
From z ~ 3 to 5, this effect is due to difference in dust redden- 
ing and age estimations, and in the case of +NEB+Lya due to 
the contribution from hya line, which can decrease the UV flux 
necessary to fit the measured fluxes. 

Relying on our previous identification of "strong" nebu- 
lar emitters and "weak" nebular emitters (Section |4.1| ), we are 
able to check the consistency of star formation rate estimation. 
Indeed, since emission lines are produced by the strong UV flux 
from OB stars in H II regions, we should find a higher SFR for 
"strong" nebular emitters than for "weak" nebular emitters. As 
shown in Figure 15 the REF model does not reproduce such a 
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Fig. 17. Same as Figure 13 for M* and SFR. The dashed line 
represents the SFR-M* relation found in Dadd i et al. (2007) at 

z~2. 



separation, since SFRs are roughly similar for both populations 
or even showing the opposite trend to what is expected. The in- 
clusion of nebular emission with a constant SF does not lead to 
a significant difference on median SFR estimations between the 
two populations. Other models without nebular lines (DEC and 
RIS) do not provide better result than constant SF since they also 
provide a trend opposite to what is expected: larger median SFR 
for "weak" nebular emitters. On the other hand, when nebular 
emission is included the two populations are naturally separated 
in terms of SFR, as shown in Fig. 15 revealing the "strong" 
nebular emitters as objecs with a strong ongoing star formation 
episode, and "weak" emitters as a more quiescent population. 
The capacity of both the declining and rising star formation his- 
tories to distinguish these populations can be easily understood, 
since in this case at least two star formation regime are possible. 
Separating the two LBG populations we find that the median 
SFR is higher by ~ 0.6 dex (up to 0.75) for the "strong" nebular 
emitters from z ~ 3 to z ~ 5 compared to the "weak" emitters, al- 
though the typical uncertainty is relatively large (~ 0.5 dex). At 
Z ~ 6, only DEC/RIS+NEB+Lya allow to retrieve the expected 
trend, showing that the presence of just one emission line like 
Lya can have a large impact on parameter estimation ( Schaerer 
|etal.|2011) . 

We now explore the SFR-M1500 relation, illustrated in Figure 
16 for the z ~ 4 sample using three different models sets. For the 
constant star formation (REF) model, the SFR-M1500 match with 



Kennicutt calibration (Kennicutt 1998), taking into account the 
effect of dust extinction. As explained in |Kennicutt| ( |1998| , the 
relation is valid for galaxies with continuous star formation over 
time scales of 10 8 years. SFR/L1500 will be significantly higher 
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Fig. 18. Same as Figure 13 for specific SFR. 



in bursty galaxies with a decreasing SF and a short timescale, or 
simply galaxies younger than 10 8 years. This will lead to signifi- 
cantly higher SFR than those given by the Kennicutt relation, re- 
gardless of dust reddening. The SFRs found cover a large range 
of possible values, strongly dependin g of both S FH and fitting or 
not with nebular emission (see Table 
only rising and declining SF with nel 



AJ||A^| and[A3|). Again, 
are able to 



ibular emission, 

separate the two populations previously identified as "strong" 
and "weak" nebular emitters, since they naturally separate these 
groups into higher and lower SFR galaxies. 



The SFR as a function M* is plotted in Figure 17 for the 
z ~ 4 sample. The figure shows that for constant star formation 
(R EF model), w e find a relation compatible to that found at z ~ 2 
by |Daddi et al. (2007 ) with a relatively small dispersion, and no 
significant difference if we consider nebular emission. For de- 
creasing SF, our results remain compatible with the relation at 
z ~ 2 but with a very large dispersion, which can be explained 
by the large range of timescale, while for rising SF, the star for- 
mation rates are systematically higher than those expected from 
the SFR-mass relation derived at z ~ 2. We note that for rising 
and decreasing star formation the galaxies seems to be separated 
in two groups, where actively star forming galaxies show higher 
SFRs than expected from the Daddi et aL] ( |2007| ) relation, and 
a group of more quiescent galaxies is compatible with this rela- 
tion. These groups correspond again largely to those previously 
identified as 'weak" nebular emitters and "strong" nebular emit- 
ters. 

The specific star formation rate (sSFR = SFR/M*) is plotted 
for z ~ 4 as function of stellar mass in Figure 18 For all models 
it decreases on average with increasing M* and with decreasing 
redshift. The relation is fairly similar among all models, but de- 
creasing and rising star formation histories lead to higher sSFR 
values, by more than 1 dex. This increase is significant compared 
to the typical errors, which range from ~ 0.2 dex for models 
with constant SF to ~ 0.6 for decreasing and rising SFHs. For 
decreasing and rising SF, the presence of Lya leads to slightly 
lower SFR and higher M+ (except at z ~ 6 where this trend 
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is reversed), which explains the lower sSFR when compared to 
models assuming no Lya emission. Comparing declining star 
formation histories to others yield lower sSFR for some galaxies. 
Both with the DEC+NEB(+Lya) and RIS+NEB(+Lyof) models, 
"strong" nebular emitters have a lower median M* and higher 
sSFR than the "weak" nebular emitters. This decrease of sSFR 
with M+ can have a physical explanation sin ce the star forma- 
tion efficiency decreases with stellar mass (eg. Kauffmann et al. 
|2003] ) . 

4.3.8. Metallicity 

Metallicity is the parameter least constrained by our SED fits. 
Indeed, for individual objects the 68% confidence interval for 
all samples covers basically the three metallicity values (0.02, 
0.2, 1 Z Q ) used here. Considering the median metallicity there is 
a trend for RIS+NEB(+Lya) and DEC+NEB(+Lya) models to 
show an increase of the metallicity with galaxy mass. However, 
the uncertainties are too large to provide firm conclusions. This 
is consistent with the well known fact that photometry is a poor 
metallicity indicator. 

4.4. Evolution of the physical properties from z~ 3 toz ~ 6 

We will now examine the evolution of the median physical prop- 
erties with redshift and discuss several implications. To allow a 
meaningful comparison avoiding in particular variations of the 
completeness limit and of the galaxy luminosity function with z, 
we do so in bins of absolute UV magnitude. To discuss the phys- 
ical parameters derived from models including nebular emission 
we choose the models with Lya flux set to zero (+NEB option) 
for z ~ 3-5, and the NEB+Lya option with maximum Lya emis- 
sion for z ~ 6. Within the options discussed in this paper, this 
choice describes best the trend of increasing strength of Lya with 
redshift as observed from spectroscopic su rveys (eg. |Ando et al. 
2006 l|Quchi et al.|2008HStark et al.|20I0| ), other studies (Blanc 
et al.|2011||Hayes et al.|2011 ), and from o ur models allowing for 
varying Lya strength ( Schaerer et al.|20lT] ). 

In the following, while we plot parameters as a function of 
M1500, in bins of magnitude, we notice that at each redshift, there 
is a very small number of object in the brightest bin, while the 
faintest bin falls beyond the completeness limit. It is therefore 
more appropriate to consider mainly the three intermediate bins 
to examine trends of the physical properties with UV magnitude. 

4.4.1. Age 
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The evolution of median ages with z is shown in Figure 
The median age increases with decreasing redshift for the four 
models, typically by one order of magnitude, with no signifi- 
cant trend with Mjjy. Although declining and rising SFs lead 
to very similar ages at the highest redshift, they start to diverge 
somewhat with decreasing redshift. We found that the age dif- 
ferences are essentially driven by the group of quiescent galax- 
ies with "weak" nebular emission. Indeed, fitting active galaxies 
with rising or decreasing SF leads to very similar ages, while 
for quiescent galaxies the ages are similar at z ~ 6 to those for 
rising SF, but start to differ at lower redshift. This is explained 
both by nebular emission and intrinsically aging population with 
decreasing redshift. Quiescent galaxies with declining star for- 
mation prefer systematically short timescales r (median at 10 
Myr at each redshift) to minimize the contribution of nebular 
emission, which is achieved with ages significantly older than r 



(at least t/r ^ 2). Indeed, for such short timescales a significant 
B aimer break can then easily be obtained with only slightly older 
ages. On the other hand, for rising SF the only way to minimize 
the contribution of nebular emission (i.e. the equivalent widths 
of lines) is by choosing relatively old ages. To fit SEDs both de- 
void of strong emission lines and with a strong B aimer break, 
even older ages are, however, needed. 
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Fig. 19. Median ages given in bin of absolute magnitude at 
1500 A (no correction for dust) from z ~ 3 to 6 for four models. 
Squares illustrate constant SFR, downward (upward) triangles 
decreasing (rising) SFHs. Black symbols stand for models with- 
out nebular emission, red symbols for models including nebu- 
lar emission (NEB), and blue symbols for NEB+Lya. The error 
bars correspond to the 68% confidence limits of the probability 
distribution in each bin. Dashed lines show the amount of time 
spanned since the previous redshift bin. 



Figure 19 leads one to conclude that for all the models, ex- 
cept the one with constant SFR and no nebular emission (REF), 
the bulk of the LBG populations at each redshift are dominated 
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by young galaxies that were not present or visible at the previ- 
ous redshift. For a constant (REF+NEB) and rising (RIS+NEB) 
star formation, these results imply that an important number of 
galaxies formed in the interval from one redshift to another, 
i.e. in ~ 300-400 Myr. For decreasing star formation, due to 
the presence of actively star-forming and quiescent galaxies, we 
can interpret this result with a scenario with episodes of active 
star formation which are followed by more quiescent episodes. 
In this scenario, age becomes a poorly constrained parameter 
since the underlying older stellar population, formed in previ- 
ous episodes of strong activity, can be dominated by a newly 
emerged population. Obviously the absolute ages also depend on 
the SF timescale. For example, fixing the timescale to r = 100 
(300) Myr for the DEC+NEB models leads to no significant dif- 
ferences on the median values of the stellar mass, SFR or red- 
dening, whereas the median age increases by a factor ~ 2 (3) for 
both active and quiescent galaxies. 



4.4.2. Reddening 
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Fig. 20. Same as Figure 19 for Ay. 
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Fig. 21. Same as Figure 
model for z ~ 3 to z 
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_ for M 150 o and A v for DEC+NEB 
5 and DEC+NEB +Ly a model for z ~ 6. 
The red line is a linear fit among the whole probability distribu- 
tion function. 



As shown in Fig. [20| the median dust attenuation decreases 
with increasing redshift for all the four models considered here. 
Inclusion of nebular emission for a constant star formation does 
not provide significant difference, while decreasing and rising 
star formation histories lead respectively to +0.1-0.2 and +0.3- 
0.4 mag in reddening compared to constant SFR. S ome studies 
(eg. |Bouwens et al.||2009| |Castellano et al.||2012] ) shown that 



there is a trend of decreasing J3 with increasing absolute UV 
magnitude, which induce a similar trend on M1500 -Ay, not seen 
in Figure 20 This trend should be found whatever is the actual 
J3-Ay relation since both classical relation ( [Bouwe ns et al. 2 009 ) 
and the relation found at z ~ 4 in this study (Section 4.3.4[ ) can 
not change the sign of the slope. 

To avoid statistical effects (bins with small number of ob- 
jects), we use the whole probability distribution at each redshift 
to determine precisely the relation between A v and M uv , with a 
result shown in Figure 21 Interestingly, the only way to obtain 
at each redshift slopes consistent with the relation between (3 and 
M uv ( [Bouwens et al. 2009 1 Castel lano et al.|2012| ), is to consider 
DEC+NEB model from z ~ 3 to z ~ 5, and DEC+NEB +Ly a at 
z ~ 6. For example, rising SF leads to positive slope at z ~ 5, as 
constant SF at z ~ 6 (both for any model considered). 

As already mentioned in Sect. |4.3.4| the median UV atten- 
uation obtained from our models with rising and declining star 
formation histories and nebular emission is higher than predicted 
from conventional methods relying on the average UV slope. For 
all models, reddening evolves similarly with redshift: it increases 
with decreasing redshift. 

Unfortunately, while photometry in the filters commonly 
used to determine J3 (Bouwe ns et al.||20TT] ) is available for the 
whole sample at z ~ 4, this is not the case at z ~ 5 and 6, where 
less than 20% of the necessary fluxes are available for the V- 
and /-drop samples. To circumvent this shortcoming we use the 
fluxes predicted by the best- fit model in these filters. At z ~ 5, 
using the equations of B ouwens et al.| ( |20lT] ) to derive J3, we also 
find a deviation from the classical A v -J3 relation, although less 
important than at z ~ 4, with an intermediate relation between 
the classical relation and that we found, while at z ~ 6, the result 
is near from equations [T]|2] and [3] but does not allow us to distin- 
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guish among them since at low Ay they are very near and there 
is only few objects at z ~ 6 with significantly high extinction. 
We have to notice that we use a / band from VLT/IS AAC while 
jBouwens et al.| ( |20TT] ) use data from Hubble/WFC3, which can 
lead to some differences. 



4.4.3. Star formation rate 

The behavior of the instantaneous SFR as a function of the UV 



magnitude and for all redshifts is shown in Figs. 22 and 23 
where we have separated the sample in the two groups, LBGs 
with "strong" and "weak" nebular emission (also referred to as 
'active" and "quiescent" galaxies here), identified earlier. For 
quiescent galaxies, all the models lead to very similar median 
SFRs, while for active galaxies, rising and declining star for- 
mation histories lead to significantly higher median SFRs, as 
already discussed. For each group we do not observe a signif- 
icant change of the SFR with redshifts. The SFR-M UV relation 
does not evolve with redshift since median age increases with 
decreasing redshift, while median reddening decreases, so the 
two effects cancel out on average. 

In the next section we discuss some implications from these 
results on the star formation history. 



4.4.4. Stellar mass and implications for star formation history 



Star k et aL] ([2009 ) suggest to explore the star formation history 
of LBGs through the evolution of the stellar mass in bins of UV 
magnitude, since the (observed/uncorrected) UV magnitude is 
the most reliable tool at high redshift to track star formation, 
and since stellar mass is the most reliable parameter, i.e. the one 
least dependent on different assumptions (SFH, metallicity, dust) 
( jFinlator et al.|2007[|Yabe et al.|2009) . Despite the fa ct that this 
study challenges this last assumption (Section 4.3.6 ), the M+- 
M1500 relation remains a useful tool able to test/falsify some 
sc enarii, mainly con stant star formation history. As explained 
Stark et aL] (2009), if the bulk of galaxies formed stars with a 



constant star formation rate over sufficiently long time we would 
expect to see a systematic increase of the normalisation of the 
M*-Mi5oo relation with cosmic time, i.e. decreasing redshift. If 
no such change is observed, the SFR cannot be constant, at least 
not over more than At ^ 300-400 Myr, the time corresponding 
to Az ~ 1 between our samples. In other words a non-evolution 
of the M+-M 1500 relation with redshift would require other star 
formation histories or the repeated emergence of new galaxies 
dominating at each redshift. 

The predicted relation between mass and UV magnitude ob- 
tained from all model sets and for all redshifts is shown in Fig. 
24| Although obviously the absolute stellar masses depend on 
the model assumptions, all models yield essentially no evolution 
of the M+-M 1500 relation between z ~ 5 to 3, but some shift be- 
tween z ~ 6 and 5. In particular from redshift 5 to 3, models 
with constant star formation (both with or without nebular emis- 
sion, i.e. REF+NEB or REF) do not yield an evolution of stellar 
mass which would be consistent with the assumption of con- 
stant SFR. Even considering the median ages obtained from the 
fits (Figure [19]), it seems difficult to reconcile the picture of con- 
stant star formation with the derived parameters. In contrast, the 
evolution of M*-Mi5oo from z ~ 6 to 5, for REF model, is fully 
compatible with what is expected from a constant star formation, 
while for the same model with nebular emission, we still do not 
see the expected evolution. However, for this latter case, median 
ages allow to assume that at each redshift, young LBGs domi- 
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Fig. 22. Same as Figure 19 for SFR and for active galaxies. The 



dashed line shows the Kennicutt relation ( Kenn Sutt|1998| ). 



nate samples. We conclude that constant SF over long timescales 
is not compatible with the data. 

Assuming a constant star formation and no strong evolution 
of the dust contain or its geometric distribution with cosmic time, 
we expect galaxies to evolve at constant M uv . With this SFH, 
the evolution of UV luminosity function should be mainly due 
to the emergence of new galaxies. For M uv = -20.5, the number 
of gala xies per M pc 3 increases by a factor ~ 3 between z ~ 6 and 
z ~ 4 fBouwens et al.|2007] ), which means that at least one third 
of LBGs seen at z ~ 4 must have LBGs at z ~ 6 as progenitors 
at this magnitude. Indeed, at z ~ 4 and for M uv ^ -20.5, our 
sample shows one third of the LBGs having age > 700 Myr, for 
REF+NEB model, which correspond to the time spans between 
z ~ 6 to z ~ 4. We predict the median stellar mass of this galax- 
ies at z ~ 4, using our parameter estimation at z ~ 6: the median 
stellar mass should be ^ 10 10-11 M*, while our sample of galax- 
ies at z ~ 4 sufficently older to be descendants of z ~ 6 galaxies 
have a median stellar mass ^ 10 9 55 . This difference of ~ 0.6 
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Fig. 23. Same as Figure 19 for SFR and for quiescent galaxies. 
The dashed line shows the Kennicutt relation. 



dex is three times larger than typical uncertainty on stellar mass 
estimation with REF+NEB model, leading us to conclude that 
despite a better consistency of the data with constant SF when 
nebular emission is taken into account, there is still significative 
discrepancies. 

What about the rising star formation history adopted here? 
It is straightforward to compute how the stellar mass and SFR 
of each galaxy evolves with time, both forward and backward. 
Starting from the best-fit for the 1000 MC simulations for each 
object, we determine how they would evolve in the future. 
Taking for e.g. LBGs at z ~ 5 which have a median stellar mass 
of 10 89 M and a median SFR ~ 30 M yr _1 , we find that they 
would have a median stellar mass of 10 1L4 M© and median SFR 
~ 2000 M Q yr" 1 (corresponding to M uv - -26 assuming the 
Kennicutt relation) after ~ 400 Myr at z ~ 4. Obviously such 
extreme objects are not seen, certainly not in large numbers! To 
hide most of them from the LBG selection a strong increase of 
dust attenuation with time would be necessary, and in this case it 
should be fairly easy to find this large populations as IR/sub-mm 



Fig. 24. Same as Figure 19 for stellar mass. In each panel, we 
plot the M*-Mi5oo found at z ~ 4 for more convenien t com - 
parison with other studies ( [Stark et al.||2009HLee et al.||2011 ). 
Black solid line: relationship for REF model, dashed black line: 
REF+NEB model, dashed red line: DEC+NEB, and dashed- 
dotted red line: RIS+NEB. 



galaxies. More likely this discrepancy implies that the average 
ris ing star formation hi story adopted here from the simulations 
of Finlator et al.| ( |2011| ) are not representative for typical z ^ 6 
galaxies, that the growth does not continue significantly beyond 
the current ages of the LBGs, or that the growth is predicted to be 
too fast, or that a combination of these arguments apply. Rising 
SFHs with varying timescales will be explored in a future pa- 
per. In any case, the study of Papo vich et al.| ( |2011| ) shows for 
the observed number counts at z ~ 3-8 to be compatible with a 
cosmologically average rising star formation history, the growth 
of the SFR has to be fairly slow. Their SFR(7) oc (t/r) a with 
a = 1.7 + 0.2 and r = 180 + 40 Myr corresponds to an in- 
crease of the SFR by a factor ~ 2 per Az = 1 from redshift 6 to 
3, a much slower growth than predicted during the first ~ 100- 



400 Myr of the star formation history of Finlator et al. (2011 ). 
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Fig. 25. Same as Figure 17 for three models (Top: REF, cen- 
tre: DEC+NEB, bottom: RIS+NEB) at z ~ 3 (Left) and z ~ 6 
(Right). For z ~ 6, models with declining and rising SF include 
hya. 



A consequence of episodic SF is that age becomes a poorly 
constrained parameter. It is also difficult to place constraints on 
the timescale of activity and inactivity, and indeed, to determine 
if estimated parameters at different redshift are consistent with 
this scenario. |Lee et al" ( |2009a| ) provide a higher limit on star 
formation activity of 350 Myr, based on observe d UV luminosity 
function and clustering at z ~ 4 and z ~ 5, and Wyithe & Loeb 
(2011) found that at high-redshift (z > 6) starburst timescale 
is set by the lifetime of massive stars, by comparing different 
assumptions on supernova feedback and observed evolution in 
galaxy size and UV luminosity function. These two results pro- 
mote short timescales, which seems to be confirmed by our re- 
sults, while uncertainties remain large. Furthermore, the result 



of |Wyithe & Loeb| (2011 ) implies decreasing timescale with in- 
creasing redshift. Figure [25] show the difference in SFR-M* re- 
lation between z ~ 3 and z ~ 6: at z ~ 6, active and quiescent 
galaxies form two clearly separate groups while at z ~ 3, the two 
groups populate all the intermediate state. We can interpret this 
result with an evolving star formation timescale, with at z ~ 6 a 
shorter timescale in comparison with z ~ 3 leading to galaxies 
evolving very rapidly from an active to a quiescent state. 

Age becomes more or less a tracer of recent star formation 
activity and do not allow to determine if observed high redshift 
LBGs are progenitors of observed low redshift LBGs. In this 
scenario, the number of free parameter remains too large to con- 
clude: duration of star-formation activity, duration of inactivity, 
the possibility for LBGs to evolve into a state with no star for- 
mation activity. Furthermore, some studies ( Bouche et al.|2010t 



Wuy ts et al. 2011 ) suggest other SFHs, like exponentially in- 
creasi ng or delayed s t ar form ation. These scenarios will be stud- 
ied in |Schaerer etaE1 ( [20T2] ). 



Hence applying the average rising SFH of |Papovich et al.|p01 1| ) 
to our individual galaxieq^] should yield results fairly close to 
those obtained from our models with constant SFR, including 
for their inability to reproduce the observed diversity of emis- 
sion line strengths traced by the (3.6-4.5)yum color at z ~ 3.8-5. 
We therefore conclude that both rapidly and slowly rising star 
formation histories over long time scales (At ^ 100-200Myr) 
are not appropriate to describe individual galaxies. Some mech- 
anism turning off star formation or leading to episodic phases 
appears to be required. 

For declining SF, or other episodic star formation histories, it 
becomes very difficult to connect galaxy populations at different 
redshifts and to draw conclusions from such a comparison. This 
is of course due to strong changes in UV luminosity with time, 
and possibly also to incorrect estimates of some parameters due 
to the "outshining" effect, where the glare of young massive stars 
can hide the pro perties of older stellar population. For example, 
|Papovich et al.| ( |2001| ) estimate that an hypothetical old stellar 
population could contain up to ^ 3 - 8 times the stellar mass of 
the young stars that dominate the observed SED. 

However, |Stark et aL] ( |2009| ) predict, under the assumption 
of an episodic star formation model, the presence of massive 
objects with low star-forming activity. Since active and quiescent 
galaxies have similar distribution in stellar mass from z ~ 5 to 
z ~ 3, and significant difference in median SFR (~0.6 dex), we 
can interpret this result as a confirmation of this latter prediction. 



2 It must be recognized that the average SFH derived by Papovich 
et al. (2011 ) corresponds to a cosmologically averaged history, which a 
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Fig. 26. Median specific star formation rate as a function of red- 
shift for four models, with 68% confidence limit based on the 
whole probability distribution function, wi th comparison with 
results from different studies at M+ = 10 95 ([Noeske et al.|2007 



Daddi et al.|2007l|Stark et al.|2009[|Reddy et al.|2012| ). Atz ~ 6, 
we show results with +NEB+Lya option for declining and rising 
SFHs. Typical errors a re ~ 0.3 dex. The da shed line shows the 
relation expected from Bouche et al. ([2010 ) for a exponentially 
increasing star formation at a fixed M* = 10 9 5 M Q . The dotted 
line, given by sSFR=2 is shown to guide the eye. 



priori does not apply to individual galaxies. 
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4.5. Specific star formation rate 

Since our study gives some elements supporting an episodic star 
formation as scenario for star formation history at high redshift, 
we compare the evolution of sSFR with results from other stud- 
ies, using the compilation from Gonzalez et al. (2010), g iven at 
fixed stellar mass M+ = 5.10 9 M G ( |Noeske e t al. 20071 |Daddi 
|e7aLl|2007l |Stark et al.|2009) . Since our different models lead 
to significative different median stellar masses (for e.g., 1O 9 6 M 
for REF model at z ~ 3 and 10 8 6 M o for DEC+NEB+Lya model 
at z ~ 6), we use for comparison median sSFR for whole sam- 
ple at each redshift. Typically other studies found no significative 
change in median stellar mass with redshift and the median value 
of sSFR at M* = 5.10 9 M© is near from the median value for the 



whole sample (e.g. |Gonzalez et al.|2010| ). Due to completeness 
limit, these values can be considered as lower limit since we ob- 
serve a trend of increasing sSFR with decreasing stellar mass for 
all models, albeit this trend is moderate with constant star forma- 
tion (see Figure [T8]). However, we can not exclude the presence 
at these redshift of high star forming galaxies enshrouded with 
dust, so our results have to be considered with appropriate cau- 



tion. Our results are shown in Figure 26 

While our results with a constant star formation seems to be 
consistent with Stark et al. (2009) results for z ~ 4 and z ~ 5, we 
find a higher median sSFR at z ~ 6. Looking at IRAC channels, 
galaxies at z ~ 6 with high sSFR (~ 20 Gyr -1 ) have on average 2 
to 3 channels with no detection (or no data), while galaxies with 
lower sSFR (median ~ 1.7 Gyr -1 ) have typically no more than 
1 channel with no detection. No detection of LBGs in optical 
(rest- frame) leads to lower stellar mass estimation (~ 1 dex), 
while estimated SFR stays similar. This explains higher sSFR 
found at z ~ 6 with constant star formation. The inclusion of 
nebular emission to constant star formation (REF+NEB) leads 
to higher sSFRs, and evolution compati ble with an expone ntially 
increasing star formation proposed by |Bouche et al. p010| ). 

Our results with declining and rising star formation including 
nebular emission differ significantly from previous studies, with 
a higher median sSFR and a possible trend of increasing sSFR 
with redshift. While studies without considering nebular emis- 
sion lead to conclude that star formation seems to be drive by 
different principles below and above z ~ 2, our results seem that 
star formation history can be more homogeneous among cosmic 
time than suggested by previous studies and models which not 
include nebular emission. 



5. Discussion 

5.1. Do we obtain realistic ages? 

The dynamical timescale is a lower limit for age estimation 
of a galaxy for reasons of causality. Dynamical timescale is 



^dyn - 2rhi/cr and using different studies ( 


Bouwens et al.|2004 


Ferguson et al.||2004[|Douglas et al.||2010 


) providing values 



and |Pettini et al.| ( |1998| ) providing velocity dispersion at z ~ 3, 
and assuming no evolution of cr, we estimate that t& yn evolves 
from ~ 40 Myr at z ~ 3 to ~ 20 Myr at z ~ 6. Dynamical 



timescale is also found to scale with (1 + z) 2 (e.g. |Wyithe & 
Loeb||2011| ), which implies a variation of a factor ~ 3 between 



redshift 3 to 6. 

Since we have not imposed a lower limit for age estima- 
tion for both declining and rising SF, both models lead to a sig- 
nificant number of galaxies with age below typical dynamical 
timescale, mainly for active galaxies. Indeed median ages for 
this latter population is almost always close to the dynamical 



limitation in age. Since age estimation is dependent from star 
formation timescale for declining SF, we test some fixed values 
of r to examine effect on age and other parameters estimation. 
With r = 100 Myr, we do not observe any significant change in 
median values of the different parameters, except for age, which 
increases by a factor ~ 2.5 for both active and quiescent galaxies, 
at each redshift. However, looking at the age distribution, there 
is still a significant number of galaxies with age < f<jyn- Imposing 
a lower limit ~ 40 Myr leads to an explanation of the previous 
result: nebular emission of active galaxies seems to be correctly 
fitted only with very young ages. Indeed, using our sample at z 
e [3.8,5], with 3.6yum and 4.5/mi fluxes measured, both declin- 
ing and rising models are not any more able to produce signifi- 
cantly better fit of 3.6//m-4.5yum color (in comparison with REF 
model), which is correlated to EW(Ha). Associated meaner val- 
ues (defined in the label of Figure [4]) reach 1.4 and 1.2, respec- 
tively for DEC+NEB and RIS+NEB models, while the quality 
of the entire SEDs remain almost unchanged. 

While this discrepancy between a significant fraction of our 
estimated ages and dynamical timescale can lead to question our 
study, we remind two elements. First, declining and rising mod- 
els lead to higher uncertainty on age, typically of a factor ~ 3 
when nebular emission is taken into account, which can allow to 
reconcile almost 80% of our samples with dynamical timescale. 
Second, there is no measured velocity dispersion of nebular lines 
at z > 3, since this measurement is confronted to the limit of cur- 
rent facilities. Furthermore, we use typi cal velocity disper sion 
measured at z ~ 3 (~ 80 km.s" 1 ) from |Pettini et al.| { |l998V but 
one object of this study (among five) shows a velocity disper- 
sion more than two times larger that this latter value, which also 
decreases t^ yn by a factor two for this object. 

We conclude that apparent discrepancies between dynami- 
cal timescale and estimated ages in this study are considerably 
reduced when uncertainties on age estimation, and observations 
limitations are taken into account. 



5.2. Comparison with other studies 

During the last decade, various papers provide analysis of the 
properties of LBGs properties at high redshift. Since we use a 
large range of star formation histories, we are able to compare 
our results with several of them, carried out for various LBG 
samples between z ~ 3 and 6. The main assumptions made for 
the SED fits in these papers, the size of their galaxy samples, 
and our corresponding model for comparison are listed in Table 
[2] To allow straightforward comparisons we will use the me- 
dian/mean values derived from our probability distribution func- 
tions. Except stated otherwise, we determine the mean over our 
entire samples. Although derived from the same pdfs, the val- 
ues from our models quoted her e do therefore not correspond to 
valu es listed in Ta bles |A.lf jA3] (which are media n values). 
|Shapley et al.| p001, hereafter S01) and |Papovich et al. 



( 2001 , hereafter P01) have studied 74 and 33 LBGs respectively, 
all at z ~ 3. The sample of S01 is restricted to brighter objects 
(Muv £ -20.7) than our z ~ 3 sample, so we compare their re- 
sults only with the results for a subsample of our objects with 
the same magnitude limit. Note that both studies include pho- 
tometry up to the K-band, but not beyond. For both studies, we 
find very similar results for the median ages (S01 obtain ~ 350, 
P01 ~ 70 Myr) and the age distribution when the correspond- 
ing models without nebular emission are used. This is not sur- 
prising since age is mainly constrained by the presence of the 
B aimer Break, which both studies can constrain from Ks and 
J (or H) band data. We also find a similar median stellar mass 
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Table 2. Summary of other studies in the literature. Col. 3 indicates their assumed star formation histories, col. 4 the extinction law, 
col. 5 the number of galaxies, and col. 6 indicates our model used for comparison. 



Redshift 


Authors a 


SFH b 


ext. law 


N 


Comparison model 


z~3 


SOI 


constant SFR 


Calzetti 


74 


REF 


z~3 


P01 


exp. declining 


Calzetti 


33 


DEC 


z~4 


PG07 


exp. declining 


Calzetti & SMC 


47 


DEC 


z~4 


Lll 


constant SFR 


Calzetti 


6 C 


REF 


z~5 


S07 


SSP/constant/exp. declining 


Calzetti 


14 


DEC 


z~5 


V07 


constant SFR 


SMC 


21 


REF 


z~5 


Y09 


constant SFR 


Calzetti 


105 


REF 


z~6 


Y06 


SSP/constant SFR 


Calzetti d 


53 


DEC 


z~6 


E07 


SSP/constant/exp. declining 


Calzetti 


17 


DEC 



a SOI: Sha pley et al.|f2 001 ), P0 1:|Papovich"eTaLlpOOT]>, PG07:|Pentericci et al.|p007)>, Lll:|Lee et al.| ( [20TT] ), S07: |Stark et aT] ( [2007l >, 
V07: |Verma et al.| ( |2007| ), Y09: |Yabe et al.| ( |2009| >'Y06: |Yan et al.| ( |2006| ), and E07? lEyles et al.| ( |2007l L 
b SSP: single stellar population, burst. 
c : stack of 1913 ob jects in 6 UV magnitude bins. 
Yan et al. (2006} assume A v = 0. For a fair comparison, we run DEC model with the same assumption. 



compared to P01 (3.10 9 M©), while we obtain mean stellar mass 
lower by ~ 0.2 dex compared to S01 (2.10 10 M Q ). Since P01 
provide best-fit parameters for their sample, we can see that sig- 
nificant discrepancies appear on dust attenuation, with for both 
studies extinction factors ~ 2 times larger in comparison with 
our results. A plausible explanation is that since longer wave- 
lengths are less sensitive to reddening, we can better constrain 
this parameter with IRAC data. Furthermore, S01 used a BC96 
population synthesis which can lead to increase reddening up to 
0.2 mag in Ay in comparizon with CB01. However, S01 also in- 
dicated that this difference in population synthesis should lead 
to significantly older ages, what we do not see. 

IPentericci etaL] ( [20071 hereafter PG07) and |Lee et aLlpOTT] 
hereafter Lll) provide analysis of 47 and 1913 LBGs at z ~ 4, 
with indidivual SED fitting in the first case, while Lll used a 
stacking procedure in UV magnitude bins. The large sample of 
Lll has a lower M uv - -21 A3, so as at z ~ 3, we take care to 
do an appropriate comparison. Our results are similar with those 
from LI 1 in stellar mass, age and dust extinction, while we have 
only 17 LBGs with M uv < -21.43. A discrepancy appears with 
the results from PG07, since we find a lower mean stellar mass 
by ~ 0.4 dex, while age, dust extinction and SFRs are similar. 

Stellar masses at z ~ 5 and 6 are also found in good agree- 
ment with previous studies, when the same/similar assumptions 
are used. Concretely, at z ~ 5, V07 found a median stellar 
mass of 2.10 9 M o (3.6.10 9 M o for our study), and Y09 and S07 
found a mean stellar mass of 4.1.10 9 and 7.9.10 9 M o respec- 
tively (while we find 9.9. 10 9 and 1.2.1O 1O M with REF and DEC 
model respectively). At z ~ 6, Y06 and E07 found 9.6. 10 9 and 
1.6.10 10 M o respectively, compared to our masses of 9.9. 10 9 and 
1.1.10 10 M o . The consistency of our results with other studies, 
taking into account typical uncertainty of ~ 0.15 dex for REF 
model and ~0.2 dex for DEC model, confirm that stellar mass is 
the most reliable parameter, since different assumptions on star 
formation history and also extinction law lead to similar results 



(cf. |Papovich et al.|200H|Verma et al.|2007||Yabe et al.|2009 ). 

Ages of z ~ 5-6 LBGs obtained in the literature agree over - 
all with our results, except for the study of Verma et al. ( |2007 ). 
Indeed, compared to their young median age of 25 Myr we find 
255 Myr, a difference which cannot be explained by uncertainty 
(~ 0.15 dex for REF model and ~ 0.4 dex for DEC model). 
A possible explanation for the discrepancy with V07 could be 
the age-reddening degeneracy. However, V07 found a median 
Ay = 0.3 mag, while we find 0.2 mag, which cannot explain the 



difference on the median age estimate. The origin of this differ- 
ence remains unclear. On the other hand, Y09 and S07 found 
median ages of 25 Myr and 288 Myr (compared to 52 Myr and 
320 Myr from our study) at z ~ 5, and Y06 and E07 obtain 
290 Myr and 400 Myr (262 Myr and 190 Myr) at z ~ 6, in 
good agreement with our results when the same assumptions are 
made. 

Other difference appears on the reddening estimates: while 
Y09 found a mean Ay ~ 0.9 mag at z ~ 5, we obtain 0.4 mag 
for the corresponding constant SFR model (typical uncertainty 
~ 0.1 mag). At z ~ 6 E07 found no reddening on average, which 
we obtain Ay = 0.4 for the DEC model (typical uncertainty ~ 
0.25 mag). Our result at z ~ 5 with the DEC model (mean A v = 
0.4) is consistent with S07, who found the same mean reddening. 
The differences on dust reddening estimation can have several 
explanations: the limited size of the samples in comparison with 
ours, or the lack of deep NIR and IRAC photometry able to put 
stronger constraints both on reddening and age. However, these 
latter arguments can not be used to explain the discrepancy with 
Y09, and we did not found a satisfactory explanation. 

These discrepancies on reddening must obviously lead to dif- 
ferences on the derived star formation rates. Indeed, at z ~ 5 
Y09 found a mean SFR of 141M yr _1 and V07 found a median 
SFR of 40 M yr\ while we find ~ 5OAf yr _1 and \5M Q yr l 
from the models corresponding best to their assumptions. The 
difference on reddening estimation with Y09 explains the dif- 
ference on SFR estimation, while the result of V07 can be ex- 
plained by their very young median age, whi ch leads to a sig - 
nificant deviation from the Kennicutt relation (Kennicutt 1998 ), 
and thus to a higher SFR. S07 do not provide SFR estimations 
but Y09 fits the parameters of S07: they found a median value 
of 2OM yr _1 , consistent with our result. However, these refitted 
parameters differ significantly from original the result of S07 for 
both stellar mass and age, which casts some doubt on the ho- 
mogeneity of this comparison. For z ~ 6, Y06 and E07 found 
mean SFR < 10M o yr _1 . While our modeling with similar as- 
sumptions as Y06 lead to a similar result (7M yr _1 ), our SFR 
estimation differs significantly from E07, with an higher SFR 
(~ 80M o yr _1 ). This latter difference can be explained by the 
higher estimated dust reddening, since E07 estimated that z ~ 6 
galaxies are mainly dust free. 

Overall, for identical assumptions, we find good agreement 
with other studies in general (exception: young ages of V07 and 
high dust extinction of Y09). But when nebular emission is in- 



21 



S. de Barros et al.: Impact of nebular emission at high redshift 



eluded, the results change quite significantly, as we have shown 
here. 



5.2.1 . Evolution of the mass - UV magnitude relation? 



Stark et al. (2009 ) and Gonzalez et al. (2011 ) cover a range of 
redshift between z ~ 4 to z ~ 6 (up to z ~ 7 for the latter). 
Both studies present the M*-Muv relation, which can directly 
be compared with our results in Figure [12] and [24] As already 
discussed in Sect. 4.3.6 we find quite similar results as Gonzalez 
et al. ( 2011| ), although our model for exponentially declining 
and for the rising star formation history suggest somewhat lower 
masses at high UV luminosities. Overall our relations remain, 
however, within the scatter of ±0.5 indicated by their work. 



Both |Stark et al.| ( [2009] ) and |Gonzalez et aT] ([20TTJ find no 
evolution of the mass - UV magnitude relation with redshift. 
In contrast, our results seem to indicate a change of the M+- 
M uv relation between z ~ 5 and 6, as shown in Fig. [27] (cf. 
also Sect. |4.4.4| ). The main reason for this change is due to our 
finding of relatively young ages for z ~ 6 galaxies, whic h implies 
a lower M/L uv ratio. This effect was also noticed by |McLure| 
et al.| ( 20TT] ) for their z ~ 6-8 sample. A preliminary reanalysis 
of their sample with the same method and assumptions used in 
the present paper appears to confirm our finding from the z ~ 6 
sample. 
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Fig. 27. Composite probability distribution of M1500 and M+ 
for REF and DEC+NEB(+LyoO models at z ~ 5 and z ~ 6. 
The black dashed line represents the M*-Mi5oo trend found by 
jGonzalez et al.| ( |2011| ) and black dotted lines shows a scatter of 
±0.5 dex. The points overlaid show the median value proper- 
ties for each object in the sample, black dots for "weak" nebular 
emitters and white dots for "strong" nebular emitters. The over- 
laid contour indicates the 68% integrated probabilities on the en- 
semble properties measured from the centroid of the distribution. 
Red dots are median values of M+ in Muv bins. 



5.2.2. LBGs with strong emission lines 

An important finding from our quantitative modeling of LBGs 
with spectral models including nebular lines is the distinction 
of two separate categories of galax ies id entified as "strong" and 
"weak" nebular emitters (cf. Sect. |4.1| ) Our work has revealed 
these categories from comparisons of the fit quality for models 
with and without nebular emission, and interestingly we found 
approximately the same fraction of objects (~ 2/3 strong versus 
1/3 weak emitters) at each redshift. 

As already mentioned above, two other studies of z ~ 4- 
5 LBGs have previously found similar objects with strong Ha 
emission identified by their excess in the 3.6 yum filter with re- 
spect to 4.5 /mi. Indeed, out of a sample of 74 LBGs with spec- 
troscopic redshifts between 3.8 and 5 Shim et al. (201 1) found at 
least 6 5% of galax i es sho wing a 3.6 /im excess attributed to Ha. 
Earlier Yabe et al. ( |2009| ) had already noted a 3.6 yum excess for 
70% of their z ~ 5 LBG sample of ~ 100 galaxies, attributed to 
the same effect. Obviously our study finds the same result and a 
very similar percentage of galaxies. Our work carries this result 
further by revealing the existence of these two LBG categories 
among all the samples studied here, ranging from U-drops to 
i-drops (i.e. z ~ 3-6), and by extending this result to fainter ob- 
jects, for which no spectroscopy is currently available. 

Besides this import agreement, our results differ, however, 
on severals points from those of Shim et al. ( |2011| ). For exam- 
ple, these authors conclude from a comparison of the estimated 
Ha equivalent width (obtained from the 3.6 fim excess) and ages 
from broad-band SED fits, that 60% of their so-called Ha emit- 
ters are forming stars at a relatively constant rate, whereas the 
rest prefer more "bursty" star formation. Our models yield nearly 
opposite results (cf. Sect. 4.2 ). Since the effects of nebular emis- 
sion affect in particular the estimated ages, the models of Shim 



et al. (2011 ), which neglect these effects in their SED modeling, 
cannot give consistent physical parameters for their strong Ha 
emitters. Similarly, whereas |Shim et aL ( |2011| ) conclude for a 
preference for the SMC extinction law, our calculations for the 
B-drop sample using this law show that the vast majority of ob- 
jects is better fit with the Calzetti attenuation law, when nebular 
lines are taken into account. Finally, from our models we see no 
need for a top-heavy IMF or extremely low metallicities, sug- 
gested as p ossible causes of the strong Ha emission by |Shim 
et al. ( |20lT ). In any case, the Ha star formation rates and the Ha 

are comparable to 



equivalent widths derived [Shim et al.| ( [201 1 
ours. 



5.3. Remaining uncertainties 

5.3.1. Uncertainties affecting individual objects 

Our study illustrates in detail the way various physical param- 
eters depend on model assumptions made for the broad band 
SED fits, and a large range of parameter space has been explored 
here. Despite these extensive investigations, the impact of some 
assumptions have not been explored in depth, and several uncer- 
tainties remain. 

The impact of different extinction laws, for example, has 
hardly been discussed here. Several papers presenting SED fits 
of LBGs at different redshifts have examined e.g. the differences 
obtained with the Calzett i attenuation law (adopted here) and th e 
SMC extinction law of ([Prevot et al.|1984[|Bouchet et al.|1985|) 

;rma| 



Among them are e.g. the work of |Papovich et al.| ( |2001| ); |Verma 
~( 2007] ); [Yabe et al.| p009 ) who study LBG samples at z 



et al. 



3 and 5. A comprehensive comparison of the impact of different 
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extinction laws (incl uding SMC, LMC, Galactic, and the Calzetti 
law) is presented in Yabe et al. (2009). To examine how differ- 
ent extinction laws modify the results from SED fitting models 
including nebular emission we have modeled the B-drop sam- 
ple with the declining star formation history (DEC+NEB model) 
adopting the SMC law of ( [Prevot et al.|19 84). Qu alitatively ou r 
results show the same trends as found by Yabe et aL] ( |2009 ), 
i.e. no big difference on stellar masses, lower Ay, older ages, 
and lower SFR. However, we find that the Calzetti attenuation 
law provides a better fit for the vast majority of galaxies. Except 
for few special cases such as the lensed galaxies cB58 and the 
Cosmic Eye st udied in detail and at different wave lengths in- 
cluding the IR ( [Baker et al.|2004[|Siana et al.|2 008), where the 
SMC law appears to be favored e.g. from the measurement of the 
IR/UV luminosity and the UV slope, it is generally though that 
the Calzetti attenuation law is applicable (at least on average) to 
high redshift star forming galaxies. 

Our models including the effects of nebular emission assume 
case B recombination to compute the strength of the hydrogen 
recombination lines (and empirical line ratios for other lines 
from He and metals) and of nebular continuum emission. In par- 
ticular we do thereby assume that all Lyman continuum pho- 
tons contribute to nebular emission, neglecting therefore possi- 
ble losses of ionizing photons due to dust inside HII regions or to 
escape of Lyman continuum photons. Furthermore, we assume 
the same extinction/attenuation for stellar and nebular emission, 
whereas observations of nearby galaxies indicate that emission 
lines usually suffer from a higher extinction than the stellar con- 
tinuum (cf. Calzetti et al.|2000|). In this sense our models maxi- 



of z ~ 2 LBGs. Along similar lines we show in Schaerer et al 



mize the effects of nebular emission, whereas in reality the lines 
could be weaker than predicted by our models. On the other hand 



we have adopted empirical line ratios compiled by Anders & 
Fritz e-v. Arvensleben| ( |2003| ) from nearby galaxies, whereas the 



conditions in distant galaxies coul d be different, l eading e.g. to 
higher excitation or stronger lines (Erb et al.|2010| ) These uncer- 
tainties are currently difficult to quantity. In any case, the broad 
band SEDs reveal quite clearly the presence of nebular lines 
and our models bracket a wide range between maximum and 
no nebular emission. Future spectroscopic observations or nar- 
row/intermediate band imaging should try to provide more de- 
tailed and accurate observational constraints on the emission line 
strengths of LBGs at high redshift. Detailed predictions from our 
models regarding individual lines will be presented elsewhere 
and can be made available on request. 

As clearly shown by our study, assumptions on the star for- 
mation history have a significant impact on the estimated phys- 
ical parameters when broad band SED fits are used. Although 
certain star formati on hi stories are found to provide better fits 
than others (cf. Sect. A2_ ), it is obvious that the simple parametri- 
sations commonly adopted in the literature and in our study can 
only be very crude representations of the true SF histories of in- 
dividual galaxies. In a companion paper (Scha erer et al.||20T2| ) 
we have explored two additional families of star formation his- 
tories, so-called delayed star formation (SFR oc texp(-t/r), and 
exponentially rising histories with variable timescales. Using ar- 
bitrary star formation histories is, however, not practical, also 
given the limited number of observational constraints. Another 
approach has e.g. been taken by Finlat oret al.| p007), who have 
used the star formation history from their cosmological hydro- 
dynamic simulations to fit observed high-z galaxies. However, 
such studies have so far been limited to a very small number of 
galaxies. Future improvements on this issue, both from obser- 



( 2012) how SFHs can be distinguished by measuring their dust 
emission with future ALMA observations and/or with measure- 
ments of emission lines. 



5.3.2. Biases and selection effects 

For obvious reasons selection effects and various biases affect 
in general studies of galaxy populations, and these need to be 
taken into account in comparisons between different samples, in 
the analysis of apparent correlations between derived physical 
parameters, and in other contexts. 

To allow meaningful comparisons with other studies we have 
presented our detailed r esult s in bins of UV magnitude (see 
data in Tables 



A.l A.2 and A3). Especially in the brightest 



vations and simulations, are certainly needed. E.g. Reddy et al. 



(2012) test different SF histories using IR and UV observations 



and faintest bins the number of galaxies is quite low, so that 
these results should be taken with care. Of course, the number 
of galaxies at different redshifts various quite strongly, affecting 
also the accuracy of the median physical properties (and confi- 
dence range) derived here. 

As in most literature studies, selection effects and biases 
have not been treated here. The impact of biases on the deter- 
mination of physical paramete rs of LBGs from b road-band SED 
fitting has e.g. been studied by |Lee etaT] ([2009b), who construct 
mock galaxy catalogs from semi-analytical models, which are 
then used to fit the simulated galaxies with a standard SED fitting 
tool. They find that stellar masses can be recovered quite well, 
whereas single-component SED fitting methods underestimate 
SFRs and overestimate ages. The differences are attributed in 
part to a "mismatch" of star formation histories between their fit- 
ting tool (which assumes exponentially declining SF) and those 
predicted by their galaxy models (which are often rising). A sim- 
ilar "template mismatch" was also identified as the main cause 
for differences in the comparison of z ~ 1.5-3 merger simula- 
tions with SED fitting results carried out by | Wuyts et al. | ( [2009] ) . 
Our models are also prone to such biases/uncertainties, and the 
role of the assumed SF histories on the derived physical pa- 
rameters has already been discussed above. However, since we 
have covered a wider range of SFHs including rising SF, our 
results may less suffer from this potential problem. The results 
from SED fits with additional SFHs (including exponentially 
declining histories with adjustable timescales) are presented in 
|SchaereretaT1 ( [20T2l ). 

Various correlations have been found between physical pa- 
rameters as the stellar mass, the star formation rate, age and oth- 
ers, both in the literature and in our study. Although our work 
emphasizes in particular the way the physical parameters depend 
on various model assumptions, i.e. a relative approach, it is im- 
portant to be aware of the selection effects which may signifi- 
cantly affect such correlations. String er et al. ( 2011| ), e.g. have 
recently examined the behavior of the SFR and the specific SFR 
with stellar mass, two important quantities discussed extensively 
in the recent literature. From their simulations of mock galaxies 
to which observational selection criteria and "standard" analy- 
sis are applied, they show how true underlying trends can be 
misrepresented. Th is study also echoes the caution expressed by 
|Dunne et"aL] ( |2009| ) on the apparent sSFR-mass relation, which 
they urge, could be severely affected by selection biases. If and to 
which extent apparent corre lations between extinction and mass, 
and age and mass (cf. Sect. 4.3.6 ) may be affected by selection 
effects and biases has not yet been addressed, to the best of our 
knowledge. A quantitative study of these effects is clearly be- 
yond the scope of the present publication. 
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6. Summary and conclusions 

We present an homogeneous study of a sample of ~ 1700 
LBGs at z ~ 3 - 6 from the GOODS-MUSIC catalog dSantmi 



|et al.|[2009| ) with deep photometry from the U band to 8 jum. 
Using a modified version of the HyperZ photometric redshift 
code which take s into account nebular emission (Schaer er~&] 
de Barros 2009), we explore a range of star formation history 
(constant, exponentially decreasing and rising). We explore the 
wide parameter space in redshift, metallicity, age, an d extinc- 
tion described by the Calzetti law (Calzetti et al .|2000| ), varying 
e-folding timescales for star formation, and whether or not neb- 
ular emission is included. The main physical parameters derived 
from our models are the stellar mass, age, reddening, star for- 
mation rate SFR and specific SFR. Furthermore our models also 
provide information on the characteristic SF timescale. 

Our method and the selected sample has been described in 
Sects. [2] and [3] The detailed model results concerning the physi- 
cal parameters, correlations among them, and describing the red- 
shift evolution of the galaxy properties have been presented in 
Sect. [4] Our main results can be summarized as follows: 

- Independently of the adopted star formation history, we find 
that at all redshifts ~ 65% of the galaxies are better fit- 
ted with nebular emission and ~ 35% without (Fig. [2]). For 
galaxies with z e [3.8,5], these two groups are clearly identi- 
fied by their 3.6yum-4.5//m color (Fig.|3]), whic h is correlated 
to strong Ha emission (cf. Shim et al. 2 011| ). Furthermore, 
this color can not be reproduced if we impose an age lim- 
itation > 50Myr (Fig. [4]). This observed color distribution 
clearly indicates the existence of galaxies with strong emis- 
sion lines and others with few or no discernible signs of 
emission lines. Our SED modeling reveals the presence of 
two LBG groups (dubbed "strong" and "weak" emitters re- 
spectively) at all redshifts studied here (z ~ 3-6), from U- 
drops to i-drops. 

- Models including the effects of nebular emission and ac- 
counting for variable (declining or rising) star formation his- 
tories naturally separate the two LBG groups according to 
current star formation rate, the group of "strong" emitters 
showing a larger SFR than the "weak" emitters (Fig. [17]). In 
a scenario of declining star formation histories these groups 
could be seen as starbursts and "post-starbursts", with age 
differences compatible with this suggestion. Models with 
constant SFR and nebular emission cannot reproduce the ob- 
served range of 3.6yum-4.5yum colors of z ~ 3.8-5 galaxies. 

- Independently of the star formation history, the inclusion of 
nebular emission leads on average to significantly younger 
ages (Fig. [5}, since nebular lines in optical (rest- frame) are 
able to mimic a B aimer break. This confirms our earlier find- 
ings (Schaerer & de Barros 2009] |2010| ) for larger samples 
and over a wider redshift range. 

- We find that the derived dust attenuation mainly depends on 
the assumed star formation history, and that the treatment of 
nebular emission does not lead to a general systematic shift. 
The largest attenuation is found with rising star formation 
histories, since these always predict very recent star forma- 
tion and hence UV em ission, as already discussed e.g. by 
( Schaerer & Pello 2005 ). In this case the inclusion of neb- 
ular emission decreases somewhat the average attenuation, 
whereas the attenuation increases for declining SFHs, and 
remains unchanged for constant SFR. 

- Based on our SED fits for 700 z ~ 4 LBGs we propose a new 
average relation between the observed UV slope J3 and the at- 
tenuation Ay. Our relation deviates somewhat from the clas- 



sical relation which assumes constant SFR and ages ^100 
Myr, and leads to a higher attenuation for a given ft slope. 
Taking into account nebular emission, the stellar masses de- 
rived from the SED fits decrease by ~ 0.4 on average, and by 
larger amounts (~ 0.4 - 0.9 dex) for LBGs from the "strong" 
emitter group (Fig. [TT]). We find a trend of increasing dust 
attenuation with stellar mass (Fig. [T4]) as already suggested 
earlier for z ~ 6-8 galaxies (Schaerer & de Barros 2010), 
and a trend of increasing age with galaxy mass (Fig. [13]). 
Given the large scatter found in the SFR-M© relation for all 
models with variable star formation histories, we also find a 
large scatter for the specific star formation rate sSFR with 
stellar mass and at all redshifts. Our favored models show a 
higher media n sSFR at z ~ 3 - 6 than derived by previous 
studies (e.g. Gonzalez et al.|201l] ) . Our results tend to indi- 
cate an increase of the median sSFR with redshift, as advo- 
cated b y several theoretical galaxy formation and evolution 
models flBouche et al.||2010HDutton et al.|20101|Weinmann| 
|etal.|2011) 

While uncertainties on r remain large, our SED fits fa- 
vor short median star formation timescales (< 300 Myr). 
Furthermore we find tentative evidence for decreasing val- 
ues of r with decreasing UV luminosity among the sample 

of "strong" emitters. 

As already shown in |Stark et al.|p009] ), constant star forma- 
tion seems to be irreconcilable with the non-evolution of the 
M+ - Mi 500 relation between z ~ 5 to 3 and with the UV 
luminosity function. 

The rising average star formation of Finlator et al.| (|2011 



used in this study cannot be representative of the typical his- 
tory of many individual LBGs over long time since the pre- 
dicted increase - if continuing into the future - is too fast. 
Indeed, this would leads to a strong increase of median stel- 
lar masses and the median SFR from one redshift to another, 
which would also imply an increase of dust attenuation to 
hide these objects from the LBG selection, since such very 
massive and strongly star-forming galaxies are not seen in 
the expected numbers. 
- The presence of two groups of LBGs identified as active and 
more quiescent galaxies, our finding of general best SED fits 
with declining star formation histories, and the consistency 
of these results with constraints on duty cy cles from cluster- 
ing studies and other theoretical arguments ( |Lee et al.|20 09a; 
Wyit he & Loeb||20TT ) all concur to consider episodic star 
formation as the scenario which best fits observations. 

Our systematic and homogeneous analysis casts new light on 
the physical properties of LBGs from z ~ 3 to 6 and p ossibly also 
to higher redshift (cf. Schaerer & de Barros 2010). Obviously, 
our results have a potentially important impact on a variety of 
questions, and several implications need to be worked out. On 
the other hand, our study also calls for new observations and 
tests. 

For example, both our preferred models (variable star forma- 
tion histories with nebular emission) imply a higher UV attenu- 
ation than what is currently derived, e.g. using the observed UV 
slope. At z ~ 4 our models typically predict an increase by a 
factor ^ 3, but smaller changes at higher redshift. Implications 
on the cosmic star formation history and related topics will be 
worked out in a separate publication. If correct, a higher UV at- 
tenuation should lead to measurable changes e.g. in the IR emis- 
sion of LBGs. A stacking analysis of the LBGs studied in this 
paper shows that the models presented here are all compatible 
with the current IR, submm, and radio observations (Greve et al. 
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2012, in preparation). Detailed predictions of the IR-mm emis- 
sion from our galaxies are presented in Schaerer et al. (2012). 
Future, more sensitive observations with ALMA should be able 
to detect individual LBGs over a wide redshift range, to deter- 
mine their attenuation, and hence also to distinguish different 
star formation histories (cf. |Shim et aL||2011[ Schaerer et al. 
2012). 

An important implication from our study is that the idea of 
simple, well-defined "star formation main sequence" with the 
majority of star- forming galaxies showing tight relation between 
SFR an d M*, suggested from other studies at lower redshift (z ~ 
0-2, cf. |Daddi et adl2007HElbaz'et al.|2007l|Noeske et al.|2007) 
may not be appropriate at high redshift (z ^ 3). Indeed a rela- 
tively small scatter is only obtained assuming constant star for- 
mation over long enough timescales 50 Myr), whereas our 
models clearly provide indications for variable star formation 
histories and episodic star formation (cf. Sects. |4.3.5| |4.2| ). A 
significant scatter is also obtained for mod els ass uming rising 
or delayed star formation histories (see Sect. |4. 3. 7] and Schaerer 



et al. 2012), which are often suggested in the recent literature. 
If the scatter found from our models in the SFR-M* relation 
decreases with decreasing redshift, converging towards the re- 
sults from other studies, remains to be seen. However, a smaller 
scatter is naturally found since constant star formation is usually 
assumed in most models or for establishing the standard SFR 
calibrations used in the literature (|Daddi et al.|2007[ |Elbaz et aL] 
50071 |Gonzalez et"aLl[20TT| |Wuyts et al,||2011> or when ana- 
lyzing stacked data, which naturally smoothes out any possible 
variation ( |Lee et aL||2010| ). Establishing more precisely the at- 
tenuation, current SFR, and star formation histories of galaxies 
is therefore crucial to shed more light on these questions. 

Related to the above mentioned scatter is also the behavior 
of the specific star formation rate (sSFR) both with galaxy mass 
and on average with redshift. Basically the same questions and 
uncertainties concerning the SFR-M* relation apply here. In any 
case it must be recognized that the sSFR is strongly dependent 
on model assumptions, and seems to show a strong dependence 
on the galaxy mass. Whereas earlier determinations of the sSFR 
at z > 3 we re considered in c onflict with recent galaxy evolu- 
tion models (Bouch ^et al.|2010 ), our results are clearly in better 
agreement with the high sSFR values and the redshift evolution 
predicted by these models. A detailed confrontation of our re- 
sults with such models and more refined ones, will hopefully 
provide further insight into galaxy formation and evolution mod- 
els at high redshift. 

Finally, it is clear that our study reveals new aspects on the 
possibly complex and variable star formation histories of high 
redshift galaxies. Whereas different arguments are found in the 
literature f avouring e.g. short duty cycles and episodic star for- 
mation (cf.|Sawicki & Yee|19981|Verma et al.|2007[|Stark et al.| 



2009} |Lee et al.||2009a[ |Wyithe & Loeb||20HT based on SED 
studies, on LBG clustering and other arguments, other studies 
favor long star formation timescales (Shapley et al. 2 00 1[ |Lee| 
et al.||2011||Shim et al.||2011| ). Our study uses for the first time 
quantitatively features probing emission lines, whose strength 
is naturally sensitive to relatively rapid variations of the recent 
SFR. That such variations have been more difficult to uncover 
before, may therefore not be surprising. Direct observations of 
the emission lines in high-z LBGs should be very useful to test 
our models and provide more stringent constraints on the impor- 
tance of nebular emission and on the star formation histories of 
distant galaxies (Schaerer et al. 2012). Other studies providing 
e.g. new measurement of LBG clustering at z > 4, searches for 
passive galaxies at high redshift, or theoretical studies on star 



formation and regulation processes can also help improving our 
understanding of these important issues closely related to key 
questions on galaxy formation and evolution. 
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Table A.l. Galaxies properties over 3 < z < 6 for constant star formation and solar metallicity set of models, in M1500 bins. For 
each parameter, we give the median value and 68% confidence limits derived from the probability distribution function. Tables will 
be available in the published version. 



Table A.l. Same as Table A.l for decreasing star formation history with variable timescale and metallicity Z. 



Table A.3. Same as Table A.l for rising star formation history with variable metallicity Z. 
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